gfx/povray/lathe.c

   1 /****************************************************************************
   2 *                   lathe.c
   3 *
   4 *  This module implements functions that manipulate lathes.
   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 * Modification from Thomas Willhalm, March 1999, used with permission
  25 *
  26 *****************************************************************************/
  27
  28 /****************************************************************************
  29 *
  30 *  Explanation:
  31 *
  32 *    The lathe primitive is defined by a set of points in 2d space which
  33 *    are interpolated by linear, quadratic, or cubic splines. The resulting
  34 *    2d curve is rotated about an axis to form the final lathe object.
  35 *
  36 *    All calculations are done in the object's (u,v,w)-coordinate system
  37 *    with the (w)-axis being the rotation axis.
  38 *
  39 *    One spline segment in the (r,w)-plane is given by the equations
  40 *
  41 *      fw(t) = Aw * t^3 + Bw * t^2 + Cw * t + Dw  and
  42 *      fr(t) = Ar * t^3 + Br * t^2 + Cr * t + Dr,
  43 *
  44 *    with the parameter t ranging from 0 to 1 and r = sqrt(u*u+v*v).
  45 *
  46 *    To intersect a ray R = P + k * D transformed into the object's
  47 *    coordinate system with the lathe object, the equations
  48 *
  49 *      (Pu + k * Du)^2 + (Pv + k * Dv)^2 = (fr(t))^2
  50 *                          (Pw + k * Dw) = fw(t)
  51 *
  52 *    have to be solved for t. For valid intersections (0 <= t <= 1)
  53 *    the corresponding k can be calculated from one of the above equations.
  54 *
  55 *    Note that the degree of the polynomal to solve is two times the
  56 *    degree of the spline segments used.
  57 *
  58 *    Note that Pu, Pv, Pw and Du, Dv, Dw denote the coordinates
  59 *    of the vectors P and D.
  60 *
  61 *  Syntax:
  62 *
  63 *    lathe
  64 *    {
  65 *      [ linear_spline | quadratic_spline | cubic_spline ]
  66 *
  67 *      number_of_points,
  68 *
  69 *      <P[0]>, <P[1], ..., <P[n-1]>
  70 *
  71 *      [ sturm ]
  72 *    }
  73 *
  74 *    Note that the P[i] are 2d vectors.
  75 *
  76 *    Linear interpolation is used by default. In this case all n Points
  77 *    are used. In the quadratic case the first point is used to determine
  78 *    the derivates at the starting point P[1]. In the cubic case
  79 *    the first and last points are used to determine the derivatives at
  80 *    the starting point P[1] and ending point P[n-2].
  81 *
  82 *    To get a closed (and smooth) curve you have make sure that
  83 *
  84 *      P[0] = P[n-1] in the linear case,
  85 *
  86 *      P[0] = P[n-2] and P[1] = P[n-1] in the quadratic case, and
  87 *
  88 *      P[0] = P[n-3] and P[1] = P[n-2] and P[2] = P[n-1] in the cubic case.
  89 *
  90 *    Note that the x coordinate of a point corresponds to the r coordinate
  91 *    and the y coordinate to the w coordinate;
  92 *
  93 *  ---
  94 *
  95 *  Ideas for the lathe were taken from:
  96 *
  97 *    P. Burger and D. Gillies, "Rapid Ray Tracing of General Surfaces
  98 *    of Revolution", New Advances in Computer Graphics, Proceedings
  99 *    of CG International '89, R. A. Earnshaw, B. Wyvill (Eds.),
 100 *    Springer, ..., pp. 523-531
 101 *
 102 *    P. Burger and D. Gillies, "Swept Surfaces", Interactive Computer
 103 *    Graphics, Addison-Wesley, 1989, pp. 376-380
 104 *
 105 *  ---
 106 *
 107 *  Jun 1994 : Creation. [DB]
 108 *
 109 *****************************************************************************/
 110
 111 #include "frame.h"
 112 #include "povray.h"
 113 #include "vector.h"
 114 #include "povproto.h"
 115 #include "bbox.h"
 116 #include "bcyl.h"
 117 #include "lathe.h"
 118 #include "polysolv.h"
 119 #include "matrices.h"
 120 #include "objects.h"
 121 #include "torus.h"
 122
 123
 124
 125 /*****************************************************************************
 126 * Local preprocessor defines
 127 ******************************************************************************/
 128
 129 /* Minimal intersection depth for a valid intersection. */
 130
 131 #define DEPTH_TOLERANCE 1.0e-4
 132
 133 /* Max. number of intersecions per spline segment. */
 134
 135 #define MAX_INTERSECTIONS_PER_SEGMENT 4
 136
 137
 138
 139 /*****************************************************************************
 140 * Local typedefs
 141 ******************************************************************************/
 142
 143
 144
 145 /*****************************************************************************
 146 * Static functions
 147 ******************************************************************************/
 148
 149 static int  intersect_lathe (RAY *Ray, LATHE *Lathe, ISTACK *Depth_Stack);
 150 static int  All_Lathe_Intersections (OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack);
 151 static int  Inside_Lathe (VECTOR point, OBJECT *Object);
 152 static void Lathe_Normal (VECTOR Result, OBJECT *Object, INTERSECTION *Inter);
 153 static LATHE *Copy_Lathe (OBJECT *Object);
 154 static void Translate_Lathe (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
 155 static void Rotate_Lathe (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
 156 static void Scale_Lathe (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
 157 static void Transform_Lathe (OBJECT *Object, TRANSFORM *Trans);
 158 static void Invert_Lathe (OBJECT *Object);
 159 static void Destroy_Lathe (OBJECT *Object);
 160
 161 static int  test_hit (LATHE *, RAY *, ISTACK *, DBL, DBL, int);
 162
 163
 164 /*****************************************************************************
 165 * Local variables
 166 ******************************************************************************/
 167
 168 static METHODS Lathe_Methods =
 169 {
 170   All_Lathe_Intersections,
 171   Inside_Lathe, Lathe_Normal,
 172   (COPY_METHOD)Copy_Lathe,
 173   Translate_Lathe, Rotate_Lathe,
 174   Scale_Lathe, Transform_Lathe, Invert_Lathe, Destroy_Lathe
 175 };
 176
 177
 178
 179 /*****************************************************************************
 180 *
 181 * FUNCTION
 182 *
 183 *   All_Lathe_Intersections
 184 *
 185 * INPUT
 186 *
 187 *   Object      - Object
 188 *   Ray         - Ray
 189 *   Depth_Stack - Intersection stack
 190 *
 191 * OUTPUT
 192 *
 193 *   Depth_Stack
 194 *
 195 * RETURNS
 196 *
 197 *   int - TRUE, if a intersection was found
 198 *
 199 * AUTHOR
 200 *
 201 *   Dieter Bayer
 202 *
 203 * DESCRIPTION
 204 *
 205 *   Determine ray/lathe intersection and clip intersection found.
 206 *
 207 * CHANGES
 208 *
 209 *   Jun 1994 : Creation.
 210 *   Oct 1996 : Changed code to include faster version. [DB]
 211 *
 212 ******************************************************************************/
 213
 214 static int All_Lathe_Intersections(OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack)
 215 {
 216   Increase_Counter(stats[Ray_Lathe_Tests]);
 217
 218   if (intersect_lathe(Ray, (LATHE *)Object, Depth_Stack))
 219   {
 220     Increase_Counter(stats[Ray_Lathe_Tests_Succeeded]);
 221
 222     return(TRUE);
 223   }
 224
 225   return(FALSE);
 226 }
 227
 228
 229
 230 /*****************************************************************************
 231 *
 232 * FUNCTION
 233 *
 234 *   intersect_lathe
 235 *
 236 * INPUT
 237 *
 238 *   Ray          - Ray
 239 *   Lathe        - Lathe
 240 *   Intersection - Lathe intersection structure
 241 *
 242 * OUTPUT
 243 *
 244 *   Intersection
 245 *
 246 * RETURNS
 247 *
 248 *   int - Number of intersections found
 249 *
 250 * AUTHOR
 251 *
 252 *   Dieter Bayer
 253 *
 254 * DESCRIPTION
 255 *
 256 *   Determine ray/lathe intersection.
 257 *
 258 *   NOTE: The curve is rotated about the y-axis!
 259 *         Order of the polynomial must not be used!
 260 *
 261 * CHANGES
 262 *
 263 *   Jun 1994 : Creation.
 264 *   Oct 1996 : Changed code to include faster version. [DB]
 265 *
 266 ******************************************************************************/
 267
 268 static int intersect_lathe(RAY *Ray, LATHE *Lathe, ISTACK *Depth_Stack)
 269 {
 270   int cnt;
 271   int found, j, n1, n2;
 272   DBL k, len, r, m, w, Dy2, r0;
 273   DBL x1[7], x2[3], y1[6], y2[2];
 274   DBL best;
 275   VECTOR P, D;
 276   BCYL_INT *intervals;
 277   LATHE_SPLINE_ENTRY *Entry;
 278
 279   /* Transform the ray into the lathe space. */
 280
 281   MInvTransPoint(P, Ray->Initial, Lathe->Trans);
 282
 283   MInvTransDirection(D, Ray->Direction, Lathe->Trans);
 284
 285   VLength(len, D);
 286
 287   VInverseScaleEq(D, len);
 288
 289   /* Test if ray misses lathe's cylindrical bound. */
 290
 291 #ifdef LATHE_EXTRA_STATS
 292   Increase_Counter(stats[Lathe_Bound_Tests]);
 293 #endif
 294
 295   if (((D[Y] >= 0.0) && (P[Y] >  Lathe->Height2)) ||
 296       ((D[Y] <= 0.0) && (P[Y] <  Lathe->Height1)) ||
 297       ((D[X] >= 0.0) && (P[X] >  Lathe->Radius2)) ||
 298       ((D[X] <= 0.0) && (P[X] < -Lathe->Radius2)))
 299   {
 300     return(FALSE);
 301   }
 302
 303   /* Get distance r0 of the ray from rotation axis (i.e. y axis). */
 304
 305   r0 = fabs(P[X] * D[Z] - P[Z] * D[X]);
 306
 307   r = D[X] * D[X] + D[Z] * D[Z];
 308
 309   if (r > 0.0)
 310   {
 311     r0 /= sqrt(r);
 312   }
 313
 314   /* Test if ray misses lathe's cylindrical bound. */
 315
 316   if (r0 > Lathe->Radius2)
 317   {
 318     return(FALSE);
 319   }
 320
 321   /* Intersect all cylindrical bounds. */
 322
 323   if ((cnt = Intersect_BCyl(Lathe->Spline->BCyl, P, D)) == 0)
 324   {
 325     return(FALSE);
 326   }
 327
 328 #ifdef LATHE_EXTRA_STATS
 329   Increase_Counter(stats[Lathe_Bound_Tests_Succeeded]);
 330 #endif
 331
 332   /* Precalculate some constants that are ray-dependant only. */
 333
 334   m = D[X] * P[X] + D[Z] * P[Z];
 335
 336   Dy2 = D[Y] * D[Y];
 337
 338   /* Step through the list of intersections. */
 339
 340   found = FALSE;
 341
 342   best = BOUND_HUGE;
 343
 344   intervals = Lathe->Spline->BCyl->intervals;
 345
 346   for (j = 0; j < cnt; j++)
 347   {
 348     /* Get current segment. */
 349
 350     Entry = &Lathe->Spline->Entry[intervals[j].n];
 351
 352     /* If we already have the best intersection we may exit. */
 353
 354     if (!(Lathe->Type & IS_CHILD_OBJECT) && (intervals[j].d[0] > best))
 355     {
 356       break;
 357     }
 358
 359     /* Init number of roots found. */
 360
 361     n1 = 0;
 362
 363     /* Intersect segment. */
 364
 365     switch (Lathe->Spline_Type)
 366     {
 367       /***********************************************************************
 368       * Linear spline.
 369       ************************************************************************/
 370
 371       case LINEAR_SPLINE:
 372
 373         /* Solve 2th order polynomial. */
 374
 375         x1[0] = Entry->C[Y] * Entry->C[Y] * r - Entry->C[X] * Entry->C[X] * Dy2;
 376
 377         x1[1] = 2.0 * (Entry->C[Y] * ((Entry->D[Y] - P[Y]) * r + D[Y] * m) - Entry->C[X] * Entry->D[X] * Dy2);
 378
 379         x1[2] = (Entry->D[Y] - P[Y]) * ((Entry->D[Y] - P[Y]) * r + 2.0 * D[Y] * m) + Dy2 * (P[X] * P[X] + P[Z] * P[Z] - Entry->D[X] * Entry->D[X]);
 380
 381         n1 = Solve_Polynomial(2, x1, y1, FALSE, 0.0);
 382
 383         break;
 384
 385       /***********************************************************************
 386       * Quadratic spline.
 387       ************************************************************************/
 388
 389       case QUADRATIC_SPLINE:
 390
 391         /* Solve 4th order polynomial. */
 392
 393         x1[0] = Entry->B[Y] * Entry->B[Y] * r - Entry->B[X] * Entry->B[X] * Dy2;
 394
 395         x1[1] = 2.0 * (Entry->B[Y] * Entry->C[Y] * r - Entry->B[X] * Entry->C[X] * Dy2);
 396
 397         x1[2] = r * (2.0 * Entry->B[Y] * (Entry->D[Y] - P[Y]) + Entry->C[Y] * Entry->C[Y]) + 2.0 * Entry->B[Y] * D[Y] * m - (2.0 * Entry->B[X] * Entry->D[X] + Entry->C[X] * Entry->C[X]) * Dy2;
 398
 399         x1[3] = 2.0 * (Entry->C[Y] * ((Entry->D[Y] - P[Y]) * r + D[Y] * m) - Entry->C[X] * Entry->D[X] * Dy2);
 400
 401         x1[4] = (Entry->D[Y] - P[Y]) * ((Entry->D[Y] - P[Y]) * r + 2.0 * D[Y] * m) + Dy2 * (P[X] * P[X] + P[Z] * P[Z] - Entry->D[X] * Entry->D[X]);
 402
 403         n1 = Solve_Polynomial(4, x1, y1, Test_Flag(Lathe, STURM_FLAG), 0.0);
 404
 405         break;
 406
 407       /***********************************************************************
 408       * Cubic spline.
 409       ************************************************************************/
 410
 411       case BEZIER_SPLINE:
 412       case CUBIC_SPLINE:
 413
 414         /* Solve 6th order polynomial. */
 415
 416         x1[0] = Entry->A[Y] * Entry->A[Y] * r - Entry->A[X] * Entry->A[X] * Dy2;
 417
 418         x1[1] = 2.0 * (Entry->A[Y] * Entry->B[Y] * r - Entry->A[X] * Entry->B[X] * Dy2);
 419
 420         x1[2] = (2.0 * Entry->A[Y] * Entry->C[Y] + Entry->B[Y] * Entry->B[Y]) * r - (2.0 * Entry->A[X] * Entry->C[X] + Entry->B[X] * Entry->B[X]) * Dy2;
 421
 422         x1[3] = 2.0 * ((Entry->A[Y] * Entry->D[Y] + Entry->B[Y] * Entry->C[Y] - Entry->A[Y] * P[Y]) * r + Entry->A[Y] * D[Y] * m - (Entry->A[X] * Entry->D[X] + Entry->B[X] * Entry->C[X]) * Dy2);
 423
 424         x1[4] = (2.0 * Entry->B[Y] * (Entry->D[Y] - P[Y]) + Entry->C[Y] * Entry->C[Y]) * r + 2.0 * Entry->B[Y] * D[Y] * m - (2.0 * Entry->B[X] * Entry->D[X] + Entry->C[X] * Entry->C[X]) * Dy2;
 425
 426         x1[5] = 2.0 * (Entry->C[Y] * ((Entry->D[Y] - P[Y]) * r + D[Y] * m) - Entry->C[X] * Entry->D[X] * Dy2);
 427
 428         x1[6] = (Entry->D[Y] - P[Y]) * ((Entry->D[Y] - P[Y]) * r + 2.0 * D[Y] * m) + Dy2 * (P[X] * P[X] + P[Z] * P[Z] - Entry->D[X] * Entry->D[X]);
 429
 430         n1 = Solve_Polynomial(6, x1, y1, Test_Flag(Lathe, STURM_FLAG), 0.0);
 431
 432         break;
 433     }
 434
 435     /* Test roots for valid intersections. */
 436
 437     while (n1--)
 438     {
 439       w = y1[n1];
 440
 441       if ((w >= 0.0) && (w <= 1.0))
 442       {
 443         if (fabs(D[Y]) > EPSILON)
 444         {
 445           k = (w * (w * (w * Entry->A[Y] + Entry->B[Y]) + Entry->C[Y]) + Entry->D[Y] - P[Y]) / D[Y];
 446
 447           if (test_hit(Lathe, Ray, Depth_Stack, k / len, w, intervals[j].n))
 448           {
 449             found = TRUE;
 450
 451             if (k < best)
 452             {
 453               best = k;
 454             }
 455           }
 456         }
 457         else
 458         {
 459           k = w * (w * (w * Entry->A[X] + Entry->B[X]) + Entry->C[X]) + Entry->D[X];
 460
 461           x2[0] = r;
 462           x2[1] = 2.0 * m;
 463
 464           x2[2] = P[X] * P[X] + P[Z] * P[Z] - k * k;
 465
 466           n2 = Solve_Polynomial(2, x2, y2, FALSE, 0.0);
 467
 468           while (n2--)
 469           {
 470             k = y2[n2];
 471
 472             if (test_hit(Lathe, Ray, Depth_Stack, k / len, w, intervals[j].n))
 473             {
 474               found = TRUE;
 475
 476               if (k < best)
 477               {
 478                 best = k;
 479               }
 480             }
 481           }
 482         }
 483       }
 484     }
 485   }
 486
 487   return(found);
 488 }
 489
 490
 491
 492 /*****************************************************************************
 493 *
 494 * FUNCTION
 495 *
 496 *   Inside_Lathe
 497 *
 498 * INPUT
 499 *
 500 *   IPoint - Intersection point
 501 *   Object - Object
 502 *
 503 * OUTPUT
 504 *
 505 * RETURNS
 506 *
 507 *   int - TRUE if inside
 508 *
 509 * AUTHOR
 510 *
 511 *   Dieter Bayer
 512 *
 513 * DESCRIPTION
 514 *
 515 *   Test if a point lies inside the lathe.
 516 *
 517 * CHANGES
 518 *
 519 *   Jun 1994 : Creation.
 520 *
 521 ******************************************************************************/
 522
 523 static int Inside_Lathe(VECTOR IPoint, OBJECT *Object)
 524 {
 525   int i, n, NC;
 526   DBL r, k, w;
 527   DBL x[4], y[3];
 528   DBL *height;
 529   VECTOR P;
 530   BCYL_ENTRY *entry;
 531   LATHE_SPLINE_ENTRY *Entry;
 532   LATHE *Lathe = (LATHE *)Object;
 533
 534   height = Lathe->Spline->BCyl->height;
 535
 536   entry = Lathe->Spline->BCyl->entry;
 537
 538   /* Transform the point into the lathe space. */
 539
 540   MInvTransPoint(P, IPoint, Lathe->Trans);
 541
 542   /* Number of crossings. */
 543
 544   NC = 0;
 545
 546   if ((P[Y] >= Lathe->Height1) && (P[Y] <= Lathe->Height2))
 547   {
 548     r = sqrt(P[X] * P[X] + P[Z] * P[Z]);
 549
 550     if ((r >= Lathe->Radius1) && (r <= Lathe->Radius2))
 551     {
 552       for (i = 0; i < Lathe->Number; i++)
 553       {
 554         Entry = &Lathe->Spline->Entry[i];
 555
 556         /* Test if we are inside the segments cylindrical bound. */
 557
 558         if ((P[Y] >= height[entry[i].h1] - EPSILON) &&
 559             (P[Y] <= height[entry[i].h2] + EPSILON))
 560         {
 561           x[0] = Entry->A[Y];
 562           x[1] = Entry->B[Y];
 563           x[2] = Entry->C[Y];
 564           x[3] = Entry->D[Y] - P[Y];
 565
 566           n = Solve_Polynomial(3, x, y, Test_Flag(Lathe, STURM_FLAG), 0.0);
 567
 568           while (n--)
 569           {
 570             w = y[n];
 571
 572             if ((w >= 0.0) && (w <= 1.0))
 573             {
 574               k  = w * (w * (w * Entry->A[X] + Entry->B[X]) + Entry->C[X]) + Entry->D[X] - r;
 575
 576               if (k >= 0.0)
 577               {
 578                 NC++;
 579               }
 580             }
 581           }
 582         }
 583       }
 584     }
 585   }
 586
 587   if (NC & 1)
 588   {
 589     return(!Test_Flag(Lathe, INVERTED_FLAG));
 590   }
 591   else
 592   {
 593     return(Test_Flag(Lathe, INVERTED_FLAG));
 594   }
 595 }
 596
 597
 598
 599 /*****************************************************************************
 600 *
 601 * FUNCTION
 602 *
 603 *   Lathe_Normal
 604 *
 605 * INPUT
 606 *
 607 *   Result - Normal vector
 608 *   Object - Object
 609 *   Inter  - Intersection found
 610 *
 611 * OUTPUT
 612 *
 613 *   Result
 614 *
 615 * RETURNS
 616 *
 617 * AUTHOR
 618 *
 619 *   Dieter Bayer
 620 *
 621 * DESCRIPTION
 622 *
 623 *   Calculate the normal of the lathe in a given point.
 624 *
 625 * CHANGES
 626 *
 627 *   Jun 1994 : Creation.
 628 *
 629 ******************************************************************************/
 630
 631 static void Lathe_Normal(VECTOR Result, OBJECT *Object, INTERSECTION *Inter)
 632 {
 633   DBL r, dx, dy;
 634   VECTOR P, N;
 635   LATHE *Lathe = (LATHE *)Object;
 636   LATHE_SPLINE_ENTRY *Entry;
 637
 638   Entry = &Lathe->Spline->Entry[Inter->i1];
 639
 640   /* Transform the point into the lathe space. */
 641
 642   MInvTransPoint(P, Inter->IPoint, Lathe->Trans);
 643
 644   /* Get distance from rotation axis. */
 645
 646   r = P[X] * P[X] + P[Z] * P[Z];
 647
 648   if (r > EPSILON)
 649   {
 650     r = sqrt(r);
 651
 652     /* Get derivatives. */
 653
 654     dx = Inter->d1 * (3.0 * Entry->A[X] * Inter->d1 + 2.0 * Entry->B[X]) + Entry->C[X];
 655     dy = Inter->d1 * (3.0 * Entry->A[Y] * Inter->d1 + 2.0 * Entry->B[Y]) + Entry->C[Y];
 656
 657     /* Get normal by rotation. */
 658
 659     N[X] =  dy * P[X];
 660     N[Y] = -dx * r;
 661     N[Z] =  dy * P[Z];
 662   }
 663   else
 664   {
 665     N[X] = N[Z] = 0.0; N[Y] = 1.0;
 666   }
 667
 668   /* Transform the normalt out of the lathe space. */
 669
 670   MTransNormal(Result, N, Lathe->Trans);
 671
 672   VNormalize(Result, Result);
 673 }
 674
 675
 676
 677 /*****************************************************************************
 678 *
 679 * FUNCTION
 680 *
 681 *   Translate_Lathe
 682 *
 683 * INPUT
 684 *
 685 *   Object - Object
 686 *   Vector - Translation vector
 687 *
 688 * OUTPUT
 689 *
 690 *   Object
 691 *
 692 * RETURNS
 693 *
 694 * AUTHOR
 695 *
 696 *   Dieter Bayer
 697 *
 698 * DESCRIPTION
 699 *
 700 *   Translate a lathe.
 701 *
 702 * CHANGES
 703 *
 704 *   Jun 1994 : Creation.
 705 *
 706 ******************************************************************************/
 707
 708 static void Translate_Lathe(OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
 709 {
 710   Transform_Lathe(Object, Trans);
 711 }
 712
 713
 714
 715 /*****************************************************************************
 716 *
 717 * FUNCTION
 718 *
 719 *   Rotate_Lathe
 720 *
 721 * INPUT
 722 *
 723 *   Object - Object
 724 *   Vector - Rotation vector
 725 *
 726 * OUTPUT
 727 *
 728 *   Object
 729 *
 730 * RETURNS
 731 *
 732 * AUTHOR
 733 *
 734 *   Dieter Bayer
 735 *
 736 * DESCRIPTION
 737 *
 738 *   Rotate a lathe.
 739 *
 740 * CHANGES
 741 *
 742 *   Jun 1994 : Creation.
 743 *
 744 ******************************************************************************/
 745
 746 static void Rotate_Lathe(OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
 747 {
 748   Transform_Lathe(Object, Trans);
 749 }
 750
 751
 752
 753 /*****************************************************************************
 754 *
 755 * FUNCTION
 756 *
 757 *   Scale_Lathe
 758 *
 759 * INPUT
 760 *
 761 *   Object - Object
 762 *   Vector - Scaling vector
 763 *
 764 * OUTPUT
 765 *
 766 *   Object
 767 *
 768 * RETURNS
 769 *
 770 * AUTHOR
 771 *
 772 *   Dieter Bayer
 773 *
 774 * DESCRIPTION
 775 *
 776 *   Scale a lathe.
 777 *
 778 * CHANGES
 779 *
 780 *   Jun 1994 : Creation.
 781 *
 782 ******************************************************************************/
 783
 784 static void Scale_Lathe(OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
 785 {
 786   Transform_Lathe(Object, Trans);
 787 }
 788
 789
 790
 791 /*****************************************************************************
 792 *
 793 * FUNCTION
 794 *
 795 *   Transform_Lathe
 796 *
 797 * INPUT
 798 *
 799 *   Object - Object
 800 *   Trans  - Transformation to apply
 801 *
 802 * OUTPUT
 803 *
 804 *   Object
 805 *
 806 * RETURNS
 807 *
 808 * AUTHOR
 809 *
 810 *   Dieter Bayer
 811 *
 812 * DESCRIPTION
 813 *
 814 *   Transform a lathe and recalculate its bounding box.
 815 *
 816 * CHANGES
 817 *
 818 *   Jun 1994 : Creation.
 819 *
 820 ******************************************************************************/
 821
 822 static void Transform_Lathe(OBJECT *Object, TRANSFORM *Trans)
 823 {
 824   Compose_Transforms(((LATHE *)Object)->Trans, Trans);
 825
 826   Compute_Lathe_BBox((LATHE *)Object);
 827 }
 828
 829
 830
 831 /*****************************************************************************
 832 *
 833 * FUNCTION
 834 *
 835 *   Invert_Lathe
 836 *
 837 * INPUT
 838 *
 839 *   Object - Object
 840 *
 841 * OUTPUT
 842 *
 843 *   Object
 844 *
 845 * RETURNS
 846 *
 847 * AUTHOR
 848 *
 849 *   Dieter Bayer
 850 *
 851 * DESCRIPTION
 852 *
 853 *   Invert a lathe.
 854 *
 855 * CHANGES
 856 *
 857 *   Jun 1994 : Creation.
 858 *
 859 ******************************************************************************/
 860
 861 static void Invert_Lathe(OBJECT *Object)
 862 {
 863   Invert_Flag(Object, INVERTED_FLAG);
 864 }
 865
 866
 867
 868 /*****************************************************************************
 869 *
 870 * FUNCTION
 871 *
 872 *   Create_Lathe
 873 *
 874 * INPUT
 875 *
 876 * OUTPUT
 877 *
 878 * RETURNS
 879 *
 880 *   LATHE * - new lathe
 881 *
 882 * AUTHOR
 883 *
 884 *   Dieter Bayer
 885 *
 886 * DESCRIPTION
 887 *
 888 *   Create a new lathe.
 889 *
 890 * CHANGES
 891 *
 892 *   Jun 1994 : Creation.
 893 *
 894 ******************************************************************************/
 895
 896 LATHE *Create_Lathe()
 897 {
 898   LATHE *New;
 899
 900   New = (LATHE *)POV_MALLOC(sizeof(LATHE), "lathe");
 901
 902   INIT_OBJECT_FIELDS(New,LATHE_OBJECT,&Lathe_Methods)
 903
 904   New->Trans = Create_Transform();
 905
 906   New->Spline_Type = LINEAR_SPLINE;
 907
 908   New->Number = 0;
 909
 910   New->Spline = NULL;
 911
 912   New->Radius1 =
 913   New->Radius2 =
 914   New->Height1 =
 915   New->Height2 = 0.0;
 916
 917   return(New);
 918 }
 919
 920
 921
 922 /*****************************************************************************
 923 *
 924 * FUNCTION
 925 *
 926 *   Copy_Lathe
 927 *
 928 * INPUT
 929 *
 930 *   Object - Object
 931 *
 932 * OUTPUT
 933 *
 934 * RETURNS
 935 *
 936 *   void * - New lathe
 937 *
 938 * AUTHOR
 939 *
 940 *   Dieter Bayer
 941 *
 942 * DESCRIPTION
 943 *
 944 *   Copy a lathe structure.
 945 *
 946 *   NOTE: The splines are not copied, only the number of references is
 947 *         counted, so that Destray_Lathe() knows if they can be destroyed.
 948 *
 949 * CHANGES
 950 *
 951 *   Jun 1994 : Creation.
 952 *
 953 *   Sep 1994 : Fixed memory leakage bug. [DB]
 954 *
 955 ******************************************************************************/
 956
 957 static LATHE *Copy_Lathe(OBJECT *Object)
 958 {
 959   LATHE *New, *Lathe = (LATHE *)Object;
 960
 961   New = Create_Lathe();
 962
 963   /* Get rid of the transformation created in Create_Lathe(). */
 964
 965   Destroy_Transform(New->Trans);
 966
 967   /* Copy lathe. */
 968
 969   *New = *Lathe;
 970
 971   New->Trans = Copy_Transform(Lathe->Trans);
 972
 973   New->Spline->References++;
 974
 975   return(New);
 976 }
 977
 978
 979
 980 /*****************************************************************************
 981 *
 982 * FUNCTION
 983 *
 984 *   Destroy_Lathe
 985 *
 986 * INPUT
 987 *
 988 *   Object - Object
 989 *
 990 * OUTPUT
 991 *
 992 *   Object
 993 *
 994 * RETURNS
 995 *
 996 * AUTHOR
 997 *
 998 *   Dieter Bayer
 999 *
1000 * DESCRIPTION
1001 *
1002 *   Destroy a lathe.
1003 *
1004 *   NOTE: The splines are destroyed if they are no longer used by any copy.
1005 *
1006 * CHANGES
1007 *
1008 *   Jun 1994 : Creation.
1009 *   Oct 1996 : Changed code to include faster version. [DB]
1010 *
1011 ******************************************************************************/
1012
1013 static void Destroy_Lathe(OBJECT *Object)
1014 {
1015   LATHE *Lathe = (LATHE *)Object;
1016
1017   Destroy_Transform(Lathe->Trans);
1018
1019   if (--(Lathe->Spline->References) == 0)
1020   {
1021     Destroy_BCyl(Lathe->Spline->BCyl);
1022
1023     POV_FREE(Lathe->Spline->Entry);
1024
1025     POV_FREE(Lathe->Spline);
1026   }
1027
1028   POV_FREE (Object);
1029 }
1030
1031
1032
1033 /*****************************************************************************
1034 *
1035 * FUNCTION
1036 *
1037 *   Compute_Lathe_BBox
1038 *
1039 * INPUT
1040 *
1041 *   Lathe - Lathe
1042 *
1043 * OUTPUT
1044 *
1045 *   Lathe
1046 *
1047 * RETURNS
1048 *
1049 * AUTHOR
1050 *
1051 *   Dieter Bayer
1052 *
1053 * DESCRIPTION
1054 *
1055 *   Calculate the bounding box of a lathe.
1056 *
1057 * CHANGES
1058 *
1059 *   Jun 1994 : Creation.
1060 *
1061 ******************************************************************************/
1062
1063 void Compute_Lathe_BBox(LATHE *Lathe)
1064 {
1065   Make_BBox(Lathe->BBox, -Lathe->Radius2, Lathe->Height1, -Lathe->Radius2,
1066     2.0 * Lathe->Radius2, Lathe->Height2 - Lathe->Height1, 2.0 * Lathe->Radius2);
1067
1068   Recompute_BBox(&Lathe->BBox, Lathe->Trans);
1069 }
1070
1071
1072
1073 /*****************************************************************************
1074 *
1075 * FUNCTION
1076 *
1077 *   Compute_Lathe
1078 *
1079 * INPUT
1080 *
1081 *   Lathe - Lathe
1082 *   P     - Points defining lathe
1083 *
1084 * OUTPUT
1085 *
1086 *   Lathe
1087 *
1088 * RETURNS
1089 *
1090 * AUTHOR
1091 *
1092 *   Dieter Bayer
1093 *
1094 * DESCRIPTION
1095 *
1096 *   Calculate the spline segments of a lathe from a set of points.
1097 *
1098 *   Note that the number of points in the lathe has to be set.
1099 *
1100 * CHANGES
1101 *
1102 *   Jun 1994 : Creation.
1103 *   Oct 1996 : Changed code to include faster version. [DB]
1104 *
1105 ******************************************************************************/
1106
1107 void Compute_Lathe(LATHE *Lathe, UV_VECT *P)
1108 {
1109   int i, i1, i2, i3, n, segment, number_of_segments;
1110   DBL x[4], xmin, xmax;
1111   DBL y[4], ymin, ymax;
1112   DBL c[3], r[2];
1113   DBL *tmp_r1;
1114   DBL *tmp_r2;
1115   DBL *tmp_h1;
1116   DBL *tmp_h2;
1117   UV_VECT A, B, C, D;
1118
1119   /* Get number of segments. */
1120
1121   switch (Lathe->Spline_Type)
1122   {
1123     case LINEAR_SPLINE:
1124
1125       number_of_segments = Lathe->Number - 1;
1126
1127       break;
1128
1129     case QUADRATIC_SPLINE:
1130
1131       number_of_segments = Lathe->Number - 2;
1132
1133       break;
1134
1135     case CUBIC_SPLINE:
1136
1137       number_of_segments = Lathe->Number - 3;
1138
1139       break;
1140
1141     case BEZIER_SPLINE:
1142
1143       number_of_segments = Lathe->Number / 4;
1144
1145       break;
1146
1147   default: /* tw */
1148
1149     number_of_segments = 0; /* tw */
1150
1151   }
1152
1153   /* Allocate segments. */
1154
1155   if (Lathe->Spline == NULL)
1156   {
1157     Lathe->Spline = (LATHE_SPLINE *)POV_MALLOC(sizeof(LATHE_SPLINE), "spline segments of lathe");
1158
1159     /* Init spline. */
1160
1161     Lathe->Spline->References = 1;
1162
1163     Lathe->Spline->Entry = (LATHE_SPLINE_ENTRY *)POV_MALLOC(number_of_segments*sizeof(LATHE_SPLINE_ENTRY), "spline segments of lathe");
1164   }
1165   else
1166   {
1167     /* This should never happen! */
1168
1169     Error("Lathe segments are already defined.\n");
1170   }
1171
1172   /* Allocate temporary lists. */
1173
1174   tmp_r1 = (DBL *)POV_MALLOC(number_of_segments * sizeof(DBL), "temp lathe data");
1175   tmp_r2 = (DBL *)POV_MALLOC(number_of_segments * sizeof(DBL), "temp lathe data");
1176   tmp_h1 = (DBL *)POV_MALLOC(number_of_segments * sizeof(DBL), "temp lathe data");
1177   tmp_h2 = (DBL *)POV_MALLOC(number_of_segments * sizeof(DBL), "temp lathe data");
1178
1179   /***************************************************************************
1180   * Calculate segments.
1181   ****************************************************************************/
1182
1183   /* We want to know the size of the overall bounding cylinder. */
1184
1185   xmax = ymax = -BOUND_HUGE;
1186   xmin = ymin = BOUND_HUGE;
1187
1188   for (i = segment = 0; segment < number_of_segments; )
1189   {
1190     i1 = i + 1;
1191     i2 = i + 2;
1192     i3 = i + 3;
1193
1194     switch (Lathe->Spline_Type)
1195     {
1196       /*************************************************************************
1197       * Linear spline (nothing more than a simple polygon).
1198       **************************************************************************/
1199
1200       case LINEAR_SPLINE:
1201
1202         /* Use linear interpolation. */
1203
1204         A[X] =  0.0;
1205         B[X] =  0.0;
1206         C[X] = -1.0 * P[i][X] + 1.0 * P[i1][X];
1207         D[X] =  1.0 * P[i][X];
1208
1209         A[Y] =  0.0;
1210         B[Y] =  0.0;
1211         C[Y] = -1.0 * P[i][Y] + 1.0 * P[i1][Y];
1212         D[Y] =  1.0 * P[i][Y];
1213
1214         /* Get maximum coordinates in current segment. */
1215
1216         x[0] = x[2] = P[i][X];
1217         x[1] = x[3] = P[i1][X];
1218
1219         y[0] = y[2] = P[i][Y];
1220         y[1] = y[3] = P[i1][Y];
1221
1222         break;
1223
1224
1225       /*************************************************************************
1226       * Quadratic spline.
1227       **************************************************************************/
1228
1229       case QUADRATIC_SPLINE:
1230
1231         /* Use quadratic interpolation. */
1232
1233         A[X] =  0.0;
1234         B[X] =  0.5 * P[i][X] - 1.0 * P[i1][X] + 0.5 * P[i2][X];
1235         C[X] = -0.5 * P[i][X]                  + 0.5 * P[i2][X];
1236         D[X] =                  1.0 * P[i1][X];
1237
1238         A[Y] =  0.0;
1239         B[Y] =  0.5 * P[i][Y] - 1.0 * P[i1][Y] + 0.5 * P[i2][Y];
1240         C[Y] = -0.5 * P[i][Y]                  + 0.5 * P[i2][Y];
1241         D[Y] =                  1.0 * P[i1][Y];
1242
1243         /* Get maximum coordinates in current segment. */
1244
1245         x[0] = x[2] = P[i1][X];
1246         x[1] = x[3] = P[i2][X];
1247
1248         y[0] = y[2] = P[i1][Y];
1249         y[1] = y[3] = P[i2][Y];
1250
1251         break;
1252
1253
1254       /*************************************************************************
1255       * Cubic spline.
1256       **************************************************************************/
1257
1258       case CUBIC_SPLINE:
1259
1260         /* Use cubic interpolation. */
1261
1262         A[X] = -0.5 * P[i][X] + 1.5 * P[i1][X] - 1.5 * P[i2][X] + 0.5 * P[i3][X];
1263         B[X] =        P[i][X] - 2.5 * P[i1][X] + 2.0 * P[i2][X] - 0.5 * P[i3][X];
1264         C[X] = -0.5 * P[i][X]                  + 0.5 * P[i2][X];
1265         D[X] =                        P[i1][X];
1266
1267         A[Y] = -0.5 * P[i][Y] + 1.5 * P[i1][Y] - 1.5 * P[i2][Y] + 0.5 * P[i3][Y];
1268         B[Y] =        P[i][Y] - 2.5 * P[i1][Y] + 2.0 * P[i2][Y] - 0.5 * P[i3][Y];
1269         C[Y] = -0.5 * P[i][Y]                  + 0.5 * P[i2][Y];
1270         D[Y] =                        P[i1][Y];
1271
1272         /* Get maximum coordinates in current segment. */
1273
1274         x[0] = x[2] = P[i1][X];
1275         x[1] = x[3] = P[i2][X];
1276
1277         y[0] = y[2] = P[i1][Y];
1278         y[1] = y[3] = P[i2][Y];
1279
1280         break;
1281
1282       /*************************************************************************
1283       * Bezier spline.
1284       **************************************************************************/
1285
1286       case BEZIER_SPLINE:
1287
1288         /* Use Bernstein interpolation. */
1289
1290         A[X] = P[i][X] - 3.0 * P[i1][X] + 3.0 * P[i2][X] -       P[i3][X];
1291         B[X] =           3.0 * P[i1][X] - 6.0 * P[i2][X] + 3.0 * P[i3][X];
1292         C[X] =           3.0 * P[i2][X] - 3.0 * P[i3][X];
1293         D[X] =                                                   P[i3][X];
1294
1295         A[Y] = P[i][Y] - 3.0 * P[i1][Y] + 3.0 * P[i2][Y] -       P[i3][Y];
1296         B[Y] =           3.0 * P[i1][Y] - 6.0 * P[i2][Y] + 3.0 * P[i3][Y];
1297         C[Y] =           3.0 * P[i2][Y] - 3.0 * P[i3][Y];
1298         D[Y] =                                                   P[i3][Y];
1299
1300         x[0] = P[i][X];
1301         x[1] = P[i1][X];
1302         x[2] = P[i2][X];
1303         x[3] = P[i3][X];
1304
1305         y[0] = P[i][Y];
1306         y[1] = P[i1][Y];
1307         y[2] = P[i2][Y];
1308         y[3] = P[i3][Y];
1309
1310         break;
1311
1312       default:
1313
1314         Error("Unknown lathe type in Compute_Lathe().\n");
1315
1316     }
1317
1318     Assign_UV_Vect(Lathe->Spline->Entry[segment].A, A);
1319     Assign_UV_Vect(Lathe->Spline->Entry[segment].B, B);
1320     Assign_UV_Vect(Lathe->Spline->Entry[segment].C, C);
1321     Assign_UV_Vect(Lathe->Spline->Entry[segment].D, D);
1322
1323     if ((Lathe->Spline_Type == QUADRATIC_SPLINE) ||
1324         (Lathe->Spline_Type == CUBIC_SPLINE))
1325     {
1326       /* Get maximum coordinates in current segment. */
1327
1328       c[0] = 3.0 * A[X];
1329       c[1] = 2.0 * B[X];
1330       c[2] = C[X];
1331
1332       n = Solve_Polynomial(2, c, r, FALSE, 0.0);
1333
1334       while (n--)
1335       {
1336         if ((r[n] >= 0.0) && (r[n] <= 1.0))
1337         {
1338           x[n] = r[n] * (r[n] * (r[n] * A[X] + B[X]) + C[X]) + D[X];
1339         }
1340       }
1341
1342       c[0] = 3.0 * A[Y];
1343       c[1] = 2.0 * B[Y];
1344       c[2] = C[Y];
1345
1346       n = Solve_Polynomial(2, c, r, FALSE, 0.0);
1347
1348       while (n--)
1349       {
1350         if ((r[n] >= 0.0) && (r[n] <= 1.0))
1351         {
1352           y[n] = r[n] * (r[n] * (r[n] * A[Y] + B[Y]) + C[Y]) + D[Y];
1353         }
1354       }
1355     }
1356
1357     /* Set current segment's bounding cylinder. */
1358
1359     tmp_r1[segment] = min(min(x[0], x[1]), min(x[2], x[3]));
1360     tmp_r2[segment] = max(max(x[0], x[1]), max(x[2], x[3]));
1361
1362     tmp_h1[segment] = min(min(y[0], y[1]), min(y[2], y[3]));
1363     tmp_h2[segment] = max(max(y[0], y[1]), max(y[2], y[3]));
1364
1365     /* Keep track of overall bounding cylinder. */
1366
1367     xmin = min(xmin, tmp_r1[segment]);
1368     xmax = max(xmax, tmp_r2[segment]);
1369
1370     ymin = min(ymin, tmp_h1[segment]);
1371     ymax = max(ymax, tmp_h2[segment]);
1372
1373 /*
1374     fprintf(stderr, "bound spline segment %d: ", i);
1375     fprintf(stderr, "r = %f - %f, h = %f - %f\n", tmp_r1[segment], tmp_r2[segment], tmp_h1[segment], tmp_h2[segment]);
1376 */
1377
1378     /* Advance to next segment. */
1379
1380     switch (Lathe->Spline_Type)
1381     {
1382       case LINEAR_SPLINE:
1383       case QUADRATIC_SPLINE:
1384       case CUBIC_SPLINE:
1385
1386         i++;
1387
1388         break;
1389
1390       case BEZIER_SPLINE:
1391
1392         i += 4;
1393
1394         break;
1395     }
1396
1397     segment++;
1398   }
1399
1400   Lathe->Number = number_of_segments;
1401
1402   /* Set overall bounding cylinder. */
1403
1404   Lathe->Radius1 = xmin;
1405   Lathe->Radius2 = xmax;
1406
1407   Lathe->Height1 = ymin;
1408   Lathe->Height2 = ymax;
1409
1410   /* Get bounding cylinder. */
1411
1412   Lathe->Spline->BCyl = Create_BCyl(Lathe->Number, tmp_r1, tmp_r2, tmp_h1, tmp_h2);
1413
1414   /* Get rid of temp. memory. */
1415
1416   POV_FREE(tmp_h2);
1417   POV_FREE(tmp_h1);
1418   POV_FREE(tmp_r2);
1419   POV_FREE(tmp_r1);
1420 }
1421
1422
1423 /*****************************************************************************
1424 *
1425 * FUNCTION
1426 *
1427 *   test_hit
1428 *
1429 * INPUT
1430 *
1431 *   Lathe       - Pointer to lathe structure
1432 *   Ray         - Current ray
1433 *   Depth_Stack - Current depth stack
1434 *   d, w, n     - Intersection depth, parameter and segment number
1435 *
1436 * OUTPUT
1437 *
1438 *   Depth_Stack
1439 *
1440 * RETURNS
1441 *
1442 * AUTHOR
1443 *
1444 *   Dieter Bayer
1445 *
1446 * DESCRIPTION
1447 *
1448 *   Test if a hit is valid and push if on the intersection depth.
1449 *
1450 * CHANGES
1451 *
1452 *   Oct 1996 : Creation.
1453 *
1454 ******************************************************************************/
1455
1456 static int test_hit(LATHE *Lathe, RAY *Ray, ISTACK *Depth_Stack, DBL d, DBL  w, int n)
1457 {
1458   VECTOR IPoint;
1459
1460   if ((d > DEPTH_TOLERANCE) && (d < Max_Distance))
1461   {
1462     VEvaluateRay(IPoint, Ray->Initial, d, Ray->Direction);
1463
1464     if (Point_In_Clip(IPoint, ((OBJECT *)Lathe)->Clip))
1465     {
1466       push_entry_i1_d1(d, IPoint, (OBJECT *)Lathe, n, w, Depth_Stack);
1467
1468       return(TRUE);
1469     }
1470   }
1471
1472   return(FALSE);
1473 }
1474
1475
1476