gfx/povray/prism.c

   1 /****************************************************************************
   2 *                   prism.c
   3 *
   4 *  This module implements functions that manipulate prisms.
   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 *    The prism primitive is defined by a set of points in 2d space which
  31 *    are interpolated by linear, quadratic, or cubic splines. The resulting
  32 *    2d curve is swept along a straight line to form the final prism object.
  33 *
  34 *    All calculations are done in the object's (u,v,w)-coordinate system
  35 *    with the (w)-axis being the sweep axis.
  36 *
  37 *    One spline segment in the (u,v)-plane is given by the equations
  38 *
  39 *      fu(t) = Au * t * t * t + Bu * t * t + Cu * t + Du  and
  40 *      fv(t) = Av * t * t * t + Bv * t * t + Cv * t + Dv,
  41 *
  42 *    with the parameter t ranging from 0 to 1.
  43 *
  44 *    To intersect a ray R = P + k * D transformed into the object's
  45 *    coordinate system with the linear swept prism object, the equation
  46 *
  47 *      Dv * fu(t) - Du * fv(t) - (Pu * Dv - Pv * Du) = 0
  48 *
  49 *    has to be solved for t. For valid intersections (0 <= t <= 1)
  50 *    the corresponding k can be calculated from equation
  51 *
  52 *      Pu + k * Du = fu(t) or Pv + k * Dv = fv(t).
  53 *
  54 *    In the case of conic sweeping the equation
  55 *
  56 *      (Pv * Dw - Pw * Dv) * fu(t) - (Pu * Dw - Pw * Du) * fv(t)
  57 *                                  + (Pu * Dv - Pv * Du) = 0
  58 *
  59 *    has to be solved for t. For valid intersections (0 <= t <= 1)
  60 *    the corresponding k can be calculated from equation
  61 *
  62 *      Pu + k * Du = (Pw + k * Dw) * fu(t) or
  63 *      Pv + k * Dv = (Pw + k * Dw) * fv(t).
  64 *
  65 *    Note that the polynomal to solve has the same degree as the
  66 *    spline segments used.
  67 *
  68 *    Note that Pu, Pv, Pw and Du, Dv, Dw denote the coordinates
  69 *    of the vectors P and D.
  70 *
  71 *  Syntax:
  72 *
  73 *    prism {
  74 *      [ linear_sweep | cubic_sweep ]
  75 *      [ linear_spline | quadratic_spline | cubic_spline ]
  76 *
  77 *      height1, height2,
  78 *      number_of_points
  79 *
  80 *      <P(0)>, <P(1)>, ..., <P(n-1)>
  81 *
  82 *      [ open ]
  83 *      [ sturm ]
  84 *    }
  85 *
  86 *    Note that the P(i) are 2d vectors.
  87 *
  88 *  ---
  89 *
  90 *  Ideas for the prism was taken from:
  91 *
  92 *    James T. Kajiya, "New techniques for ray tracing procedurally
  93 *    defined objects", Computer Graphics, 17(3), July 1983, pp. 91-102
  94 *
  95 *    Kirk ...
  96 *
  97 *  ---
  98 *
  99 *  May 1994 : Creation. [DB]
 100 *
 101 *****************************************************************************/
 102
 103 #include "frame.h"
 104 #include "povray.h"
 105 #include "vector.h"
 106 #include "povproto.h"
 107 #include "bbox.h"
 108 #include "matrices.h"
 109 #include "objects.h"
 110 #include "polysolv.h"
 111 #include "prism.h"
 112
 113
 114
 115 /*****************************************************************************
 116 * Local preprocessor defines
 117 ******************************************************************************/
 118
 119 /* Minimal intersection depth for a valid intersection. */
 120
 121 #define DEPTH_TOLERANCE 1.0e-4
 122
 123 /* |Coefficients| < COEFF_LIMIT are regarded to be 0. */
 124
 125 /*#define COEFF_LIMIT 1.0e-20 */
 126
 127 #define COEFF_LIMIT 1.0e-16 /*changed by CEY 11/18/95 */
 128
 129 /* Part of the prism hit. */
 130
 131 #define BASE_HIT   1
 132 #define CAP_HIT    2
 133 #define SPLINE_HIT 3
 134
 135
 136
 137 /*****************************************************************************
 138 * Local typedefs
 139 ******************************************************************************/
 140
 141
 142
 143 /*****************************************************************************
 144 * Static functions
 145 ******************************************************************************/
 146
 147 static int intersect_prism (RAY *Ray, PRISM *Prism, PRISM_INT *Intersection);
 148 static int in_curve (PRISM *Prism, DBL u, DBL v);
 149 static int test_rectangle (VECTOR P, VECTOR D, DBL x1, DBL y1, DBL x2, DBL y2);
 150 static int   All_Prism_Intersections (OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack);
 151 static int   Inside_Prism (VECTOR point, OBJECT *Object);
 152 static void  Prism_Normal (VECTOR Result, OBJECT *Object, INTERSECTION *Inter);
 153 static PRISM *Copy_Prism (OBJECT *Object);
 154 static void  Translate_Prism (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
 155 static void  Rotate_Prism (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
 156 static void  Scale_Prism (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
 157 static void  Transform_Prism (OBJECT *Object, TRANSFORM *Trans);
 158 static void  Invert_Prism (OBJECT *Object);
 159 static void  Destroy_Prism (OBJECT *Object);
 160
 161
 162 /*****************************************************************************
 163 * Local variables
 164 ******************************************************************************/
 165
 166 static METHODS Prism_Methods =
 167 {
 168   All_Prism_Intersections,
 169   Inside_Prism, Prism_Normal,
 170   (COPY_METHOD)Copy_Prism,
 171   Translate_Prism, Rotate_Prism,
 172   Scale_Prism, Transform_Prism, Invert_Prism, Destroy_Prism
 173 };
 174
 175
 176
 177 /*****************************************************************************
 178 *
 179 * FUNCTION
 180 *
 181 *   All_Prism_Intersections
 182 *
 183 * INPUT
 184 *
 185 *   Object      - Object
 186 *   Ray         - Ray
 187 *   Depth_Stack - Intersection stack
 188 *
 189 * OUTPUT
 190 *
 191 *   Depth_Stack
 192 *
 193 * RETURNS
 194 *
 195 *   int - TRUE, if a intersection was found
 196 *
 197 * AUTHOR
 198 *
 199 *   Dieter Bayer
 200 *
 201 * DESCRIPTION
 202 *
 203 *   Determine ray/prism intersection and clip intersection found.
 204 *
 205 * CHANGES
 206 *
 207 *   May 1994 : Creation.
 208 *
 209 ******************************************************************************/
 210
 211 static int All_Prism_Intersections(OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack)
 212 {
 213   int i, max_i, Found;
 214   PRISM_INT *Inter;
 215   VECTOR IPoint;
 216
 217   Found = FALSE;
 218
 219   Inter = ((PRISM *)Object)->Intersections;
 220
 221   max_i = intersect_prism(Ray, (PRISM *)Object, Inter);
 222
 223   if (max_i)
 224   {
 225     for (i = 0; i < max_i; i++)
 226     {
 227       if ((Inter[i].d > DEPTH_TOLERANCE) && (Inter[i].d < Max_Distance))
 228       {
 229         VEvaluateRay(IPoint, Ray->Initial, Inter[i].d, Ray->Direction);
 230
 231         if (Point_In_Clip(IPoint, Object->Clip))
 232         {
 233           push_entry_i1_i2_d1(Inter[i].d,IPoint,Object,Inter[i].t,Inter[i].n,Inter[i].w,Depth_Stack);
 234
 235           Found = TRUE;
 236         }
 237       }
 238     }
 239   }
 240
 241   return(Found);
 242 }
 243
 244
 245
 246 /*****************************************************************************
 247 *
 248 * FUNCTION
 249 *
 250 *   intersect_prism
 251 *
 252 * INPUT
 253 *
 254 *   Ray   - Ray
 255 *   Prism - Prism
 256 *   Int   - Prism intersection structure
 257 *
 258 * OUTPUT
 259 *
 260 *   Int
 261 *
 262 * RETURNS
 263 *
 264 *   int - Number of intersections found
 265 *
 266 * AUTHOR
 267 *
 268 *   Dieter Bayer
 269 *
 270 * DESCRIPTION
 271 *
 272 *   Determine ray/prism intersection.
 273 *
 274 *   NOTE: Order reduction cannot be used.
 275 *
 276 * CHANGES
 277 *
 278 *   May 1994 : Creation.
 279 *
 280 ******************************************************************************/
 281
 282 static int intersect_prism(RAY *Ray, PRISM *Prism, PRISM_INT *Intersection)
 283 {
 284   int i, j, n;
 285   DBL k, u, v, w, h, len;
 286   DBL x[4], y[3];
 287   DBL k1, k2, k3;
 288   VECTOR P, D;
 289   PRISM_SPLINE_ENTRY Entry;
 290
 291   /* Don't test degenerate prisms. */
 292
 293   if (Test_Flag(Prism, DEGENERATE_FLAG))
 294   {
 295     return(FALSE);
 296   }
 297
 298   Increase_Counter(stats[Ray_Prism_Tests]);
 299
 300   /* Init intersection structure. */
 301
 302   Intersection->d =
 303   Intersection->w = 0.0;
 304   Intersection->n =
 305   Intersection->t = 0;
 306
 307   /* Transform the ray into the prism space */
 308
 309   MInvTransPoint(P, Ray->Initial, Prism->Trans);
 310
 311   MInvTransDirection(D, Ray->Direction, Prism->Trans);
 312
 313   VLength(len, D);
 314
 315   VInverseScaleEq(D, len);
 316
 317   /* Test overall bounding rectangle. */
 318
 319 #ifdef PRISM_EXTRA_STATS
 320   Increase_Counter(stats[Prism_Bound_Tests]);
 321 #endif
 322
 323   if (((D[X] >= 0.0) && (P[X] > Prism->x2)) ||
 324       ((D[X] <= 0.0) && (P[X] < Prism->x1)) ||
 325       ((D[Z] >= 0.0) && (P[Z] > Prism->y2)) ||
 326       ((D[Z] <= 0.0) && (P[Z] < Prism->y1)))
 327   {
 328     return(FALSE);
 329   }
 330
 331 #ifdef PRISM_EXTRA_STATS
 332   Increase_Counter(stats[Prism_Bound_Tests_Succeeded]);
 333 #endif
 334
 335   /* Number of intersections found. */
 336
 337   i = 0;
 338
 339   /* What kind of sweep is used? */
 340
 341   switch (Prism->Sweep_Type)
 342   {
 343     /*************************************************************************
 344     * Linear sweep.
 345     **************************************************************************/
 346
 347     case LINEAR_SWEEP :
 348
 349       if (fabs(D[Y]) < EPSILON)
 350       {
 351         if ((P[Y] < Prism->Height1) || (P[Y] > Prism->Height2))
 352         {
 353           return(FALSE);
 354         }
 355       }
 356       else
 357       {
 358         if (Test_Flag(Prism, CLOSED_FLAG))
 359         {
 360           /* Intersect ray with the cap-plane. */
 361
 362           k = (Prism->Height2 - P[Y]) / D[Y];
 363
 364           if ((k > DEPTH_TOLERANCE) && (k < Max_Distance))
 365           {
 366             u = P[X] + k * D[X];
 367             v = P[Z] + k * D[Z];
 368
 369             if (in_curve(Prism, u, v))
 370             {
 371               Intersection[i].t   = CAP_HIT;
 372               Intersection[i++].d = k / len;
 373             }
 374           }
 375
 376           /* Intersect ray with the base-plane. */
 377
 378           k = (Prism->Height1 - P[Y]) / D[Y];
 379
 380           if ((k > DEPTH_TOLERANCE) && (k < Max_Distance))
 381           {
 382             u = P[X] + k * D[X];
 383             v = P[Z] + k * D[Z];
 384
 385             if (in_curve(Prism, u, v))
 386             {
 387               Intersection[i].t   = BASE_HIT;
 388               Intersection[i++].d = k / len;
 389             }
 390           }
 391         }
 392       }
 393
 394       /* Intersect ray with all spline segments. */
 395
 396       if ((fabs(D[X]) > EPSILON) || (fabs(D[Z]) > EPSILON))
 397       {
 398         for (j = 0; j < Prism->Number; j++)
 399         {
 400           Entry = Prism->Spline->Entry[j];
 401
 402           /* Test spline's bounding rectangle (modified Cohen-Sutherland). */
 403
 404 #ifdef PRISM_EXTRA_STATS
 405           Increase_Counter(stats[Prism_Bound_Tests]);
 406 #endif
 407
 408           if (((D[X] >= 0.0) && (P[X] > Entry.x2)) ||
 409               ((D[X] <= 0.0) && (P[X] < Entry.x1)) ||
 410               ((D[Z] >= 0.0) && (P[Z] > Entry.y2)) ||
 411               ((D[Z] <= 0.0) && (P[Z] < Entry.y1)))
 412           {
 413             continue;
 414           }
 415
 416           /* Number of roots found. */
 417
 418           n = 0;
 419
 420           switch (Prism->Spline_Type)
 421           {
 422             case LINEAR_SPLINE :
 423
 424 #ifdef PRISM_EXTRA_STATS
 425               Increase_Counter(stats[Prism_Bound_Tests_Succeeded]);
 426 #endif
 427
 428               /* Solve linear equation. */
 429
 430               x[0] = Entry.C[X] * D[Z] - Entry.C[Y] * D[X];
 431
 432               x[1] = D[Z] * (Entry.D[X] - P[X]) - D[X] * (Entry.D[Y] - P[Z]);
 433
 434               if (fabs(x[0]) > EPSILON)
 435               {
 436                 y[n++] = -x[1] / x[0];
 437               }
 438
 439               break;
 440
 441             case QUADRATIC_SPLINE :
 442
 443 #ifdef PRISM_EXTRA_STATS
 444               Increase_Counter(stats[Prism_Bound_Tests_Succeeded]);
 445 #endif
 446
 447               /* Solve quadratic equation. */
 448
 449               x[0] = Entry.B[X] * D[Z] - Entry.B[Y] * D[X];
 450
 451               x[1] = Entry.C[X] * D[Z] - Entry.C[Y] * D[X];
 452
 453               x[2] = D[Z] * (Entry.D[X] - P[X]) - D[X] * (Entry.D[Y] - P[Z]);
 454
 455               n = Solve_Polynomial(2, x, y, FALSE, 0.0);
 456
 457               break;
 458
 459             case CUBIC_SPLINE :
 460             case BEZIER_SPLINE :
 461
 462               if (test_rectangle(P, D, Entry.x1, Entry.y1, Entry.x2, Entry.y2))
 463               {
 464 #ifdef PRISM_EXTRA_STATS
 465                 Increase_Counter(stats[Prism_Bound_Tests_Succeeded]);
 466 #endif
 467
 468                 /* Solve cubic equation. */
 469
 470                 x[0] = Entry.A[X] * D[Z] - Entry.A[Y] * D[X];
 471
 472                 x[1] = Entry.B[X] * D[Z] - Entry.B[Y] * D[X];
 473
 474                 x[2] = Entry.C[X] * D[Z] - Entry.C[Y] * D[X];
 475
 476                 x[3] = D[Z] * (Entry.D[X] - P[X]) - D[X] * (Entry.D[Y] - P[Z]);
 477
 478                 n = Solve_Polynomial(3, x, y, Test_Flag(Prism, STURM_FLAG), 0.0);
 479               }
 480
 481               break;
 482           }
 483
 484           /* Test roots for valid intersections. */
 485
 486           while (n--)
 487           {
 488             w = y[n];
 489
 490             if ((w >= 0.0) && (w <= 1.0))
 491             {
 492               if (fabs(D[X]) > EPSILON)
 493               {
 494                 k = (w * (w * (w * Entry.A[X] + Entry.B[X]) + Entry.C[X]) + Entry.D[X] - P[X]) / D[X];
 495               }
 496               else
 497               {
 498                 k = (w * (w * (w * Entry.A[Y] + Entry.B[Y]) + Entry.C[Y]) + Entry.D[Y] - P[Z]) / D[Z];
 499               }
 500
 501               /* Verify that intersection height is valid. */
 502
 503               h = P[Y] + k * D[Y];
 504
 505               if ((h >= Prism->Height1) && (h <= Prism->Height2))
 506               {
 507                 Intersection[i].t   = SPLINE_HIT;
 508                 Intersection[i].n   = j;
 509                 Intersection[i].w   = w;
 510                 Intersection[i++].d = k / len;
 511               }
 512             }
 513           }
 514         }
 515       }
 516
 517       break;
 518
 519     /*************************************************************************
 520     * Conic sweep.
 521     **************************************************************************/
 522
 523     case CONIC_SWEEP :
 524
 525       if (fabs(D[Y]) < EPSILON)
 526       {
 527         if ((P[Y] < Prism->Height1) || (P[Y] > Prism->Height2))
 528         {
 529           return(FALSE);
 530         }
 531       }
 532       else
 533       {
 534         if (Test_Flag(Prism, CLOSED_FLAG))
 535         {
 536           /* Intersect ray with the cap-plane. */
 537
 538           if (fabs(Prism->Height2) > EPSILON)
 539           {
 540             k = (Prism->Height2 - P[Y]) / D[Y];
 541
 542             if ((k > DEPTH_TOLERANCE) && (k < Max_Distance))
 543             {
 544               u = (P[X] + k * D[X]) / Prism->Height2;
 545               v = (P[Z] + k * D[Z]) / Prism->Height2;
 546
 547               if (in_curve(Prism, u, v))
 548               {
 549                 Intersection[i].t   = CAP_HIT;
 550                 Intersection[i++].d = k / len;
 551               }
 552             }
 553           }
 554
 555           /* Intersect ray with the base-plane. */
 556
 557           if (fabs(Prism->Height1) > EPSILON)
 558           {
 559             k = (Prism->Height1 - P[Y]) / D[Y];
 560
 561             if ((k > DEPTH_TOLERANCE) && (k < Max_Distance))
 562             {
 563               u = (P[X] + k * D[X]) / Prism->Height1;
 564               v = (P[Z] + k * D[Z]) / Prism->Height1;
 565
 566               if (in_curve(Prism, u, v))
 567               {
 568                 Intersection[i].t   = BASE_HIT;
 569                 Intersection[i++].d = k / len;
 570               }
 571             }
 572           }
 573         }
 574       }
 575
 576       /* Precompute ray-only dependant constants. */
 577
 578       k1 = P[Z] * D[Y] - P[Y] * D[Z];
 579
 580       k2 = P[Y] * D[X] - P[X] * D[Y];
 581
 582       k3 = P[X] * D[Z] - P[Z] * D[X];
 583
 584       /* Intersect ray with the spline segments. */
 585
 586       if ((fabs(D[X]) > EPSILON) || (fabs(D[Z]) > EPSILON))
 587       {
 588         for (j = 0; j < Prism->Number; j++)
 589         {
 590           Entry = Prism->Spline->Entry[j];
 591
 592           /* Test spline's bounding rectangle (modified Cohen-Sutherland). */
 593
 594           if (((D[X] >= 0.0) && (P[X] > Entry.x2)) ||
 595               ((D[X] <= 0.0) && (P[X] < Entry.x1)) ||
 596               ((D[Z] >= 0.0) && (P[Z] > Entry.y2)) ||
 597               ((D[Z] <= 0.0) && (P[Z] < Entry.y1)))
 598           {
 599             continue;
 600           }
 601
 602           /* Number of roots found. */
 603
 604           n = 0;
 605
 606           switch (Prism->Spline_Type)
 607           {
 608             case LINEAR_SPLINE :
 609
 610               /* Solve linear equation. */
 611
 612               x[0] = Entry.C[X] * k1 + Entry.C[Y] * k2;
 613
 614               x[1] = Entry.D[X] * k1 + Entry.D[Y] * k2 + k3;
 615
 616               if (fabs(x[0]) > EPSILON)
 617               {
 618                 y[n++] = -x[1] / x[0];
 619               }
 620
 621               break;
 622
 623             case QUADRATIC_SPLINE :
 624
 625               /* Solve quadratic equation. */
 626
 627               x[0] = Entry.B[X] * k1 + Entry.B[Y] * k2;
 628
 629               x[1] = Entry.C[X] * k1 + Entry.C[Y] * k2;
 630
 631               x[2] = Entry.D[X] * k1 + Entry.D[Y] * k2 + k3;
 632
 633               n = Solve_Polynomial(2, x, y, FALSE, 0.0);
 634
 635               break;
 636
 637             case CUBIC_SPLINE :
 638             case BEZIER_SPLINE :
 639
 640               /* Solve cubic equation. */
 641
 642               x[0] = Entry.A[X] * k1 + Entry.A[Y] * k2;
 643
 644               x[1] = Entry.B[X] * k1 + Entry.B[Y] * k2;
 645
 646               x[2] = Entry.C[X] * k1 + Entry.C[Y] * k2;
 647
 648               x[3] = Entry.D[X] * k1 + Entry.D[Y] * k2 + k3;
 649
 650               n = Solve_Polynomial(3, x, y, Test_Flag(Prism, STURM_FLAG), 0.0);
 651
 652               break;
 653           }
 654
 655           /* Test roots for valid intersections. */
 656
 657           while (n--)
 658           {
 659             w = y[n];
 660
 661             if ((w >= 0.0) && (w <= 1.0))
 662             {
 663               k = w * (w * (w * Entry.A[X] + Entry.B[X]) + Entry.C[X]) + Entry.D[X];
 664
 665               h = D[X] - k * D[Y];
 666
 667               if (fabs(h) > EPSILON)
 668               {
 669                 k = (k * P[Y] - P[X]) / h;
 670               }
 671               else
 672               {
 673                 k = w * (w * (w * Entry.A[Y] + Entry.B[Y]) + Entry.C[Y]) + Entry.D[Y];
 674
 675                 h = D[Z] - k * D[Y];
 676
 677                 if (fabs(h) > EPSILON)
 678                 {
 679                   k = (k * P[Y] - P[Z]) / h;
 680                 }
 681                 else
 682                 {
 683                   /* This should never happen! */
 684                   continue;
 685                 }
 686               }
 687
 688               /* Verify that intersection height is valid. */
 689
 690               h = P[Y] + k * D[Y];
 691
 692               if ((h >= Prism->Height1) && (h <= Prism->Height2))
 693               {
 694                 Intersection[i].t   = SPLINE_HIT;
 695                 Intersection[i].n   = j;
 696                 Intersection[i].w   = w;
 697                 Intersection[i++].d = k / len;
 698               }
 699             }
 700           }
 701         }
 702       }
 703
 704       break;
 705
 706       default:
 707
 708         Error("Unknown sweep type in intersect_prism().\n");
 709   }
 710
 711   if (i)
 712   {
 713     Increase_Counter(stats[Ray_Prism_Tests_Succeeded]);
 714   }
 715
 716   return(i);
 717 }
 718
 719
 720
 721 /*****************************************************************************
 722 *
 723 * FUNCTION
 724 *
 725 *   Inside_Prism
 726 *
 727 * INPUT
 728 *
 729 *   IPoint - Intersection point
 730 *   Object - Object
 731 *
 732 * OUTPUT
 733 *
 734 * RETURNS
 735 *
 736 *   int - TRUE if inside
 737 *
 738 * AUTHOR
 739 *
 740 *   Dieter Bayer
 741 *
 742 * DESCRIPTION
 743 *
 744 *   Test if point lies inside a prism.
 745 *
 746 * CHANGES
 747 *
 748 *   May 1994 : Creation.
 749 *
 750 ******************************************************************************/
 751
 752 static int Inside_Prism(VECTOR IPoint, OBJECT *Object)
 753 {
 754   VECTOR P;
 755   PRISM *Prism = (PRISM *)Object;
 756
 757   /* Transform the point into the prism space. */
 758
 759   MInvTransPoint(P, IPoint, Prism->Trans);
 760
 761   if ((P[Y] >= Prism->Height1) && (P[Y] < Prism->Height2))
 762   {
 763     if (Prism->Sweep_Type == CONIC_SWEEP)
 764     {
 765       /* Scale x and z coordinate. */
 766
 767       if (fabs(P[Y]) > EPSILON)
 768       {
 769         P[X] /= P[Y];
 770         P[Z] /= P[Y];
 771       }
 772       else
 773       {
 774         P[X] = P[Z] = HUGE_VAL;
 775       }
 776     }
 777
 778     if (in_curve(Prism, P[X], P[Z]))
 779     {
 780       return(!Test_Flag(Prism, INVERTED_FLAG));
 781     }
 782   }
 783
 784   return(Test_Flag(Prism, INVERTED_FLAG));
 785 }
 786
 787
 788
 789 /*****************************************************************************
 790 *
 791 * FUNCTION
 792 *
 793 *   Prism_Normal
 794 *
 795 * INPUT
 796 *
 797 *   Result - Normal vector
 798 *   Object - Object
 799 *   Inter  - Intersection found
 800 *
 801 * OUTPUT
 802 *
 803 *   Result
 804 *
 805 * RETURNS
 806 *
 807 * AUTHOR
 808 *
 809 *   Dieter Bayer
 810 *
 811 * DESCRIPTION
 812 *
 813 *   Calculate the normal of the prism in a given point.
 814 *
 815 * CHANGES
 816 *
 817 *   May 1994 : Creation.
 818 *
 819 *   Jul 1997 : Fixed bug as reported by Darko Rozic. [DB]
 820 *
 821 ******************************************************************************/
 822
 823 static void Prism_Normal(VECTOR Result, OBJECT *Object, INTERSECTION *Inter)
 824 {
 825   VECTOR P;
 826   PRISM_SPLINE_ENTRY Entry;
 827   PRISM *Prism = (PRISM *)Object;
 828   VECTOR N;
 829
 830   Make_Vector(N, 0.0, 1.0, 0.0);
 831
 832   if (Inter->i1 == SPLINE_HIT)
 833   {
 834     Entry = Prism->Spline->Entry[Inter->i2];
 835
 836     switch (Prism->Sweep_Type)
 837     {
 838       case LINEAR_SWEEP:
 839
 840         N[X] =   Inter->d1 * (3.0 * Entry.A[Y] * Inter->d1 + 2.0 * Entry.B[Y]) + Entry.C[Y];
 841         N[Y] =   0.0;
 842         N[Z] = -(Inter->d1 * (3.0 * Entry.A[X] * Inter->d1 + 2.0 * Entry.B[X]) + Entry.C[X]);
 843
 844         break;
 845
 846       case CONIC_SWEEP:
 847
 848         /* Transform the point into the prism space. */
 849
 850         MInvTransPoint(P, Inter->IPoint, Prism->Trans);
 851
 852         if (fabs(P[Y]) > EPSILON)
 853         {
 854           N[X] =   Inter->d1 * (3.0 * Entry.A[Y] * Inter->d1 + 2.0 * Entry.B[Y]) + Entry.C[Y];
 855           N[Z] = -(Inter->d1 * (3.0 * Entry.A[X] * Inter->d1 + 2.0 * Entry.B[X]) + Entry.C[X]);
 856           N[Y] = -(P[X] * N[X] + P[Z] * N[Z]) / P[Y];
 857         }
 858
 859         break;
 860
 861       default:
 862
 863         Error("Unknown sweep type in Prism_Normal().\n");
 864     }
 865   }
 866
 867   /* Transform the normalt out of the prism space. */
 868
 869   MTransNormal(Result, N, Prism->Trans);
 870
 871   VNormalize(Result, Result);
 872 }
 873
 874
 875
 876 /*****************************************************************************
 877 *
 878 * FUNCTION
 879 *
 880 *   Translate_Prism
 881 *
 882 * INPUT
 883 *
 884 *   Object - Object
 885 *   Vector - Translation vector
 886 *
 887 * OUTPUT
 888 *
 889 *   Object
 890 *
 891 * RETURNS
 892 *
 893 * AUTHOR
 894 *
 895 *   Dieter Bayer
 896 *
 897 * DESCRIPTION
 898 *
 899 *   Translate a prism.
 900 *
 901 * CHANGES
 902 *
 903 *   May 1994 : Creation.
 904 *
 905 ******************************************************************************/
 906
 907 static void Translate_Prism(OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
 908 {
 909   Transform_Prism(Object, Trans);
 910 }
 911
 912
 913
 914 /*****************************************************************************
 915 *
 916 * FUNCTION
 917 *
 918 *   Rotate_Prism
 919 *
 920 * INPUT
 921 *
 922 *   Object - Object
 923 *   Vector - Rotation vector
 924 *
 925 * OUTPUT
 926 *
 927 *   Object
 928 *
 929 * RETURNS
 930 *
 931 * AUTHOR
 932 *
 933 *   Dieter Bayer
 934 *
 935 * DESCRIPTION
 936 *
 937 *   Rotate a prism.
 938 *
 939 * CHANGES
 940 *
 941 *   May 1994 : Creation.
 942 *
 943 ******************************************************************************/
 944
 945 static void Rotate_Prism(OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
 946 {
 947   Transform_Prism(Object, Trans);
 948 }
 949
 950
 951
 952 /*****************************************************************************
 953 *
 954 * FUNCTION
 955 *
 956 *   Scale_Prism
 957 *
 958 * INPUT
 959 *
 960 *   Object - Object
 961 *   Vector - Scaling vector
 962 *
 963 * OUTPUT
 964 *
 965 *   Object
 966 *
 967 * RETURNS
 968 *
 969 * AUTHOR
 970 *
 971 *   Dieter Bayer
 972 *
 973 * DESCRIPTION
 974 *
 975 *   Scale a prism.
 976 *
 977 * CHANGES
 978 *
 979 *   May 1994 : Creation.
 980 *
 981 ******************************************************************************/
 982
 983 static void Scale_Prism(OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
 984 {
 985   Transform_Prism(Object, Trans);
 986 }
 987
 988
 989
 990 /*****************************************************************************
 991 *
 992 * FUNCTION
 993 *
 994 *   Transform_Prism
 995 *
 996 * INPUT
 997 *
 998 *   Object - Object
 999 *   Trans  - Transformation to apply
1000 *
1001 * OUTPUT
1002 *
1003 *   Object
1004 *
1005 * RETURNS
1006 *
1007 * AUTHOR
1008 *
1009 *   Dieter Bayer
1010 *
1011 * DESCRIPTION
1012 *
1013 *   Transform a prism and recalculate its bounding box.
1014 *
1015 * CHANGES
1016 *
1017 *   May 1994 : Creation.
1018 *
1019 ******************************************************************************/
1020
1021 static void Transform_Prism(OBJECT *Object, TRANSFORM *Trans)
1022 {
1023   Compose_Transforms(((PRISM *)Object)->Trans, Trans);
1024
1025   Compute_Prism_BBox((PRISM *)Object);
1026 }
1027
1028
1029
1030 /*****************************************************************************
1031 *
1032 * FUNCTION
1033 *
1034 *   Invert_Prism
1035 *
1036 * INPUT
1037 *
1038 *   Object - Object
1039 *
1040 * OUTPUT
1041 *
1042 *   Object
1043 *
1044 * RETURNS
1045 *
1046 * AUTHOR
1047 *
1048 *   Dieter Bayer
1049 *
1050 * DESCRIPTION
1051 *
1052 *   Invert a prism.
1053 *
1054 * CHANGES
1055 *
1056 *   May 1994 : Creation.
1057 *
1058 ******************************************************************************/
1059
1060 static void Invert_Prism(OBJECT *Object)
1061 {
1062   Invert_Flag(Object, INVERTED_FLAG);
1063 }
1064
1065
1066
1067 /*****************************************************************************
1068 *
1069 * FUNCTION
1070 *
1071 *   Create_Prism
1072 *
1073 * INPUT
1074 *
1075 * OUTPUT
1076 *
1077 * RETURNS
1078 *
1079 *   PRISM * - new prism
1080 *
1081 * AUTHOR
1082 *
1083 *   Dieter Bayer
1084 *
1085 * DESCRIPTION
1086 *
1087 *   Create a new prism.
1088 *
1089 * CHANGES
1090 *
1091 *   May 1994 : Creation.
1092 *
1093 ******************************************************************************/
1094
1095 PRISM *Create_Prism()
1096 {
1097   PRISM *New;
1098
1099   New = (PRISM *)POV_MALLOC(sizeof(PRISM), "prism");
1100
1101   INIT_OBJECT_FIELDS(New,PRISM_OBJECT,&Prism_Methods)
1102
1103   New->Trans = Create_Transform();
1104
1105   New->x1      =
1106   New->x2      =
1107   New->y1      =
1108   New->y2      =
1109   New->u1      =
1110   New->u2      =
1111   New->v1      =
1112   New->v2      =
1113   New->Height1 =
1114   New->Height2 = 0.0;
1115
1116   New->Number = 0;
1117
1118   New->Spline_Type = LINEAR_SPLINE;
1119   New->Sweep_Type  = LINEAR_SWEEP;
1120
1121   New->Spline = NULL;
1122
1123   New->Intersections = NULL;
1124
1125   Set_Flag(New, CLOSED_FLAG);
1126
1127   return(New);
1128 }
1129
1130
1131
1132 /*****************************************************************************
1133 *
1134 * FUNCTION
1135 *
1136 *   Copy_Prism
1137 *
1138 * INPUT
1139 *
1140 *   Object - Object
1141 *
1142 * OUTPUT
1143 *
1144 * RETURNS
1145 *
1146 *   void * - New prism
1147 *
1148 * AUTHOR
1149 *
1150 *   Dieter Bayer
1151 *
1152 * DESCRIPTION
1153 *
1154 *   Copy a prism structure.
1155 *
1156 *   NOTE: The splines are not copied, only the number of references is
1157 *         counted, so that Destray_Prism() knows if they can be destroyed.
1158 *
1159 * CHANGES
1160 *
1161 *   May 1994 : Creation.
1162 *
1163 *   Sep 1994 : fixed memory leakage [DB]
1164 *
1165 ******************************************************************************/
1166
1167 static PRISM *Copy_Prism(OBJECT *Object)
1168 {
1169   PRISM *New, *Prism = (PRISM *)Object;
1170
1171   New = Create_Prism();
1172
1173   /* Get rid of the transformation created in Create_Prism(). */
1174
1175   Destroy_Transform(New->Trans);
1176
1177   /* Copy prism. */
1178
1179   *New = *Prism;
1180
1181   New->Trans = Copy_Transform(((PRISM *)Object)->Trans);
1182
1183   New->Spline->References++;
1184
1185   Prism->Intersections = (PRISM_INT *)POV_MALLOC((Prism->Number+2)*sizeof(PRISM_INT), "prism intersection list");
1186
1187   return(New);
1188 }
1189
1190
1191
1192 /*****************************************************************************
1193 *
1194 * FUNCTION
1195 *
1196 *   Destroy_Prism
1197 *
1198 * INPUT
1199 *
1200 *   Object - Object
1201 *
1202 * OUTPUT
1203 *
1204 *   Object
1205 *
1206 * RETURNS
1207 *
1208 * AUTHOR
1209 *
1210 *   Dieter Bayer
1211 *
1212 * DESCRIPTION
1213 *
1214 *   Destroy a prism.
1215 *
1216 *   NOTE: The splines are destroyed if they are no longer used by any copy.
1217 *
1218 * CHANGES
1219 *
1220 *   May 1994 : Creation.
1221 *
1222 ******************************************************************************/
1223
1224 static void Destroy_Prism (OBJECT *Object)
1225 {
1226   PRISM *Prism = (PRISM *)Object;
1227
1228   Destroy_Transform(Prism->Trans);
1229
1230   POV_FREE(Prism->Intersections);
1231
1232   if (--(Prism->Spline->References) == 0)
1233   {
1234     POV_FREE(Prism->Spline->Entry);
1235
1236     POV_FREE(Prism->Spline);
1237   }
1238
1239   POV_FREE(Object);
1240 }
1241
1242
1243
1244 /*****************************************************************************
1245 *
1246 * FUNCTION
1247 *
1248 *   Compute_Prism_BBox
1249 *
1250 * INPUT
1251 *
1252 *   Prism - Prism
1253 *
1254 * OUTPUT
1255 *
1256 *   Prism
1257 *
1258 * RETURNS
1259 *
1260 * AUTHOR
1261 *
1262 *   Dieter Bayer
1263 *
1264 * DESCRIPTION
1265 *
1266 *   Calculate the bounding box of a prism.
1267 *
1268 * CHANGES
1269 *
1270 *   May 1994 : Creation.
1271 *
1272 ******************************************************************************/
1273
1274 void Compute_Prism_BBox(PRISM *Prism)
1275 {
1276   Make_BBox(Prism->BBox, Prism->x1, Prism->Height1, Prism->y1,
1277     Prism->x2 - Prism->x1, Prism->Height2 - Prism->Height1, Prism->y2 - Prism->y1);
1278
1279   Recompute_BBox(&Prism->BBox, Prism->Trans);
1280 }
1281
1282
1283
1284 /*****************************************************************************
1285 *
1286 * FUNCTION
1287 *
1288 *   in_curve
1289 *
1290 * INPUT
1291 *
1292 *   Prism - Prism to test
1293 *   u, v  - Coordinates
1294 *
1295 * OUTPUT
1296 *
1297 * RETURNS
1298 *
1299 *   int - TRUE if inside
1300 *
1301 * AUTHOR
1302 *
1303 *   Dieter Bayer, June 1994
1304 *
1305 * DESCRIPTION
1306 *
1307 *   Test if a 2d point lies inside a prism's spline curve.
1308 *
1309 *   We travel from the current point in positive u direction
1310 *   and count the number of crossings with the curve.
1311 *
1312 * CHANGES
1313 *
1314 *   May 1994 : Creation.
1315 *
1316 ******************************************************************************/
1317
1318 static int in_curve(PRISM *Prism, DBL u, DBL  v)
1319 {
1320   int i, n, NC;
1321   DBL k, w;
1322   DBL x[4], y[3];
1323   PRISM_SPLINE_ENTRY Entry;
1324
1325   NC = 0;
1326
1327   /* First test overall bounding rectangle. */
1328
1329   if ((u >= Prism->u1) && (u <= Prism->u2) &&
1330       (v >= Prism->v1) && (v <= Prism->v2))
1331   {
1332     for (i = 0; i < Prism->Number; i++)
1333     {
1334       Entry = Prism->Spline->Entry[i];
1335
1336       /* Test if current segment can be hit. */
1337
1338       if ((v >= Entry.v1) && (v <= Entry.v2) && (u <= Entry.u2))
1339       {
1340         x[0] = Entry.A[Y];
1341         x[1] = Entry.B[Y];
1342         x[2] = Entry.C[Y];
1343         x[3] = Entry.D[Y] - v;
1344
1345         n = Solve_Polynomial(3, x, y, Test_Flag(Prism, STURM_FLAG), 0.0);
1346
1347         while (n--)
1348         {
1349           w = y[n];
1350
1351           if ((w >= 0.0) && (w <= 1.0))
1352           {
1353             k  = w * (w * (w * Entry.A[X] + Entry.B[X]) + Entry.C[X]) + Entry.D[X] - u;
1354
1355             if (k >= 0.0)
1356             {
1357               NC++;
1358             }
1359           }
1360         }
1361       }
1362     }
1363   }
1364
1365   return(NC & 1);
1366 }
1367
1368
1369
1370 /*****************************************************************************
1371 *
1372 * FUNCTION
1373 *
1374 *   test_rectangle
1375 *
1376 * INPUT
1377 *
1378 * OUTPUT
1379 *
1380 * RETURNS
1381 *
1382 * AUTHOR
1383 *
1384 *   Dieter Bayer, July 1994
1385 *
1386 * DESCRIPTION
1387 *
1388 *   Test if the 2d ray (= P + t * D) intersects a rectangle.
1389 *
1390 * CHANGES
1391 *
1392 *   May 1994 : Creation.
1393 *
1394 ******************************************************************************/
1395
1396 static int test_rectangle(VECTOR P, VECTOR  D, DBL x1, DBL  z1, DBL  x2, DBL  z2)
1397 {
1398   DBL dmin, dmax, tmin, tmax;
1399
1400   if (fabs(D[X]) > EPSILON)
1401   {
1402     if (D[X] > 0.0)
1403     {
1404       dmin = (x1 - P[X]) / D[X];
1405       dmax = (x2 - P[X]) / D[X];
1406
1407       if (dmax < EPSILON)
1408       {
1409         return(FALSE);
1410       }
1411     }
1412     else
1413     {
1414       dmax = (x1 - P[X]) / D[X];
1415
1416       if (dmax < EPSILON)
1417       {
1418         return(FALSE);
1419       }
1420
1421       dmin = (x2 - P[X]) / D[X];
1422     }
1423
1424     if (dmin > dmax)
1425     {
1426       return(FALSE);
1427     }
1428   }
1429   else
1430   {
1431     if ((P[X] < x1) || (P[X] > x2))
1432     {
1433       return(FALSE);
1434     }
1435     else
1436     {
1437       dmin = -BOUND_HUGE;
1438       dmax =  BOUND_HUGE;
1439     }
1440   }
1441
1442   if (fabs(D[Z]) > EPSILON)
1443   {
1444     if (D[Z] > 0.0)
1445     {
1446       tmin = (z1 - P[Z]) / D[Z];
1447       tmax = (z2 - P[Z]) / D[Z];
1448     }
1449     else
1450     {
1451       tmax = (z1 - P[Z]) / D[Z];
1452       tmin = (z2 - P[Z]) / D[Z];
1453     }
1454
1455     if (tmax < dmax)
1456     {
1457       if (tmax < EPSILON)
1458       {
1459         return(FALSE);
1460       }
1461
1462       if (tmin > dmin)
1463       {
1464         if (tmin > tmax)
1465         {
1466           return(FALSE);
1467         }
1468
1469         dmin = tmin;
1470       }
1471       else
1472       {
1473         if (dmin > tmax)
1474         {
1475           return(FALSE);
1476         }
1477       }
1478
1479       /*dmax = tmax; */ /*not needed CEY[1/95]*/
1480     }
1481     else
1482     {
1483       if (tmin > dmin)
1484       {
1485         if (tmin > dmax)
1486         {
1487           return(FALSE);
1488         }
1489
1490         /* dmin = tmin; */  /*not needed CEY[1/95]*/
1491       }
1492     }
1493   }
1494   else
1495   {
1496     if ((P[Z] < z1) || (P[Z] > z2))
1497     {
1498       return(FALSE);
1499     }
1500   }
1501
1502   return(TRUE);
1503 }
1504
1505
1506
1507 /*****************************************************************************
1508 *
1509 * FUNCTION
1510 *
1511 *   Compute_Prism
1512 *
1513 * INPUT
1514 *
1515 *   Prism - Prism
1516 *   P     - Points defining prism
1517 *
1518 * OUTPUT
1519 *
1520 *   Prism
1521 *
1522 * RETURNS
1523 *
1524 * AUTHOR
1525 *
1526 *   Dieter Bayer, June 1994
1527 *
1528 * DESCRIPTION
1529 *
1530 *   Calculate the spline segments of a prism from a set of points.
1531 *
1532 * CHANGES
1533 *
1534 *   May 1994 : Creation.
1535 *
1536 ******************************************************************************/
1537
1538 void Compute_Prism(PRISM *Prism, UV_VECT *P)
1539 {
1540   int i, n, number_of_splines;
1541   int i1, i2, i3;
1542
1543   DBL x[4], xmin, xmax;
1544   DBL y[4], ymin, ymax;
1545   DBL c[3], r[2];
1546
1547   UV_VECT A, B, C, D, First;
1548
1549   /* Allocate Object->Number segments. */
1550
1551   if (Prism->Spline == NULL)
1552   {
1553     Prism->Spline = (PRISM_SPLINE *)POV_MALLOC(sizeof(PRISM_SPLINE), "spline segments of prism");
1554
1555     Prism->Spline->References = 1;
1556
1557     Prism->Spline->Entry = (PRISM_SPLINE_ENTRY *)POV_MALLOC(Prism->Number*sizeof(PRISM_SPLINE_ENTRY), "spline segments of prism");
1558   }
1559   else
1560   {
1561     /* This should never happen! */
1562
1563     Error("Prism segments are already defined.\n");
1564   }
1565
1566   /* Allocate prism intersections list. */
1567
1568   Prism->Intersections = (PRISM_INT *)POV_MALLOC((Prism->Number+2)*sizeof(PRISM_INT), "prism intersection list");
1569
1570   /***************************************************************************
1571   * Calculate segments.
1572   ****************************************************************************/
1573
1574   /* We want to know the size of the overall bounding rectangle. */
1575
1576   xmax = ymax = -BOUND_HUGE;
1577   xmin = ymin =  BOUND_HUGE;
1578
1579   /* Get first segment point. */
1580
1581   switch (Prism->Spline_Type)
1582   {
1583     case LINEAR_SPLINE:
1584
1585       Assign_UV_Vect(First, P[0]);
1586
1587       break;
1588
1589     case QUADRATIC_SPLINE:
1590     case CUBIC_SPLINE:
1591
1592       Assign_UV_Vect(First, P[1]);
1593
1594       break;
1595   }
1596
1597   for (i = number_of_splines = 0; i < Prism->Number-1; i++)
1598   {
1599     /* Get indices of previous and next two points. */
1600
1601     i1 = i + 1;
1602     i2 = i + 2;
1603     i3 = i + 3;
1604
1605     switch (Prism->Spline_Type)
1606     {
1607       /*************************************************************************
1608       * Linear spline (nothing more than a simple polygon).
1609       **************************************************************************/
1610
1611       case LINEAR_SPLINE :
1612
1613         if (i1 >= Prism->Number)
1614         {
1615           Error("Too few points in prism. Prism not closed? Control points missing?\n");
1616         }
1617
1618         /* Use linear interpolation. */
1619
1620         A[X] =  0.0;
1621         B[X] =  0.0;
1622         C[X] = -1.0 * P[i][X] + 1.0 * P[i1][X];
1623         D[X] =  1.0 * P[i][X];
1624
1625         A[Y] =  0.0;
1626         B[Y] =  0.0;
1627         C[Y] = -1.0 * P[i][Y] + 1.0 * P[i1][Y];
1628         D[Y] =  1.0 * P[i][Y];
1629
1630         x[0] = x[2] = P[i][X];
1631         x[1] = x[3] = P[i1][X];
1632
1633         y[0] = y[2] = P[i][Y];
1634         y[1] = y[3] = P[i1][Y];
1635
1636         break;
1637
1638       /*************************************************************************
1639       * Quadratic spline.
1640       **************************************************************************/
1641
1642       case QUADRATIC_SPLINE :
1643
1644         if (i2 >= Prism->Number)
1645         {
1646           Error("Too few points in prism.\n");
1647         }
1648
1649         /* Use quadratic interpolation. */
1650
1651         A[X] =  0.0;
1652         B[X] =  0.5 * P[i][X] - 1.0 * P[i1][X] + 0.5 * P[i2][X];
1653         C[X] = -0.5 * P[i][X]                  + 0.5 * P[i2][X];
1654         D[X] =                  1.0 * P[i1][X];
1655
1656         A[Y] =  0.0;
1657         B[Y] =  0.5 * P[i][Y] - 1.0 * P[i1][Y] + 0.5 * P[i2][Y];
1658         C[Y] = -0.5 * P[i][Y]                  + 0.5 * P[i2][Y];
1659         D[Y] =                  1.0 * P[i1][Y];
1660
1661         x[0] = x[2] = P[i1][X];
1662         x[1] = x[3] = P[i2][X];
1663
1664         y[0] = y[2] = P[i1][Y];
1665         y[1] = y[3] = P[i2][Y];
1666
1667         break;
1668
1669       /*************************************************************************
1670       * Cubic spline.
1671       **************************************************************************/
1672
1673       case CUBIC_SPLINE :
1674
1675         if (i3 >= Prism->Number)
1676         {
1677           Error("Too few points in prism.\n");
1678         }
1679
1680         /* Use cubic interpolation. */
1681
1682         A[X] = -0.5 * P[i][X] + 1.5 * P[i1][X] - 1.5 * P[i2][X] + 0.5 * P[i3][X];
1683         B[X] =        P[i][X] - 2.5 * P[i1][X] + 2.0 * P[i2][X] - 0.5 * P[i3][X];
1684         C[X] = -0.5 * P[i][X]                  + 0.5 * P[i2][X];
1685         D[X] =                        P[i1][X];
1686
1687         A[Y] = -0.5 * P[i][Y] + 1.5 * P[i1][Y] - 1.5 * P[i2][Y] + 0.5 * P[i3][Y];
1688         B[Y] =        P[i][Y] - 2.5 * P[i1][Y] + 2.0 * P[i2][Y] - 0.5 * P[i3][Y];
1689         C[Y] = -0.5 * P[i][Y]                  + 0.5 * P[i2][Y];
1690         D[Y] =                        P[i1][Y];
1691
1692         x[0] = x[2] = P[i1][X];
1693         x[1] = x[3] = P[i2][X];
1694
1695         y[0] = y[2] = P[i1][Y];
1696         y[1] = y[3] = P[i2][Y];
1697
1698         break;
1699
1700       /*************************************************************************
1701       * Bezier spline.
1702       **************************************************************************/
1703
1704       case BEZIER_SPLINE :
1705
1706         if (i3 >= Prism->Number)
1707         {
1708           Error("Too few points in prism. Prism not closed? Control points missing?\n");
1709         }
1710
1711         /* Use Bernstein blending function interpolation. */
1712
1713         A[X] = P[i][X] - 3.0 * P[i1][X] + 3.0 * P[i2][X] -       P[i3][X];
1714         B[X] =           3.0 * P[i1][X] - 6.0 * P[i2][X] + 3.0 * P[i3][X];
1715         C[X] =           3.0 * P[i2][X] - 3.0 * P[i3][X];
1716         D[X] =                                                   P[i3][X];
1717
1718         A[Y] = P[i][Y] - 3.0 * P[i1][Y] + 3.0 * P[i2][Y] -       P[i3][Y];
1719         B[Y] =           3.0 * P[i1][Y] - 6.0 * P[i2][Y] + 3.0 * P[i3][Y];
1720         C[Y] =           3.0 * P[i2][Y] - 3.0 * P[i3][Y];
1721         D[Y] =                                                   P[i3][Y];
1722
1723         x[0] = P[i][X];
1724         x[1] = P[i1][X];
1725         x[2] = P[i2][X];
1726         x[3] = P[i3][X];
1727
1728         y[0] = P[i][Y];
1729         y[1] = P[i1][Y];
1730         y[2] = P[i2][Y];
1731         y[3] = P[i3][Y];
1732
1733         break;
1734
1735       default:
1736
1737         Error("Unknown spline type in Compute_Prism().\n");
1738     }
1739
1740     Assign_UV_Vect(Prism->Spline->Entry[number_of_splines].A, A);
1741     Assign_UV_Vect(Prism->Spline->Entry[number_of_splines].B, B);
1742     Assign_UV_Vect(Prism->Spline->Entry[number_of_splines].C, C);
1743     Assign_UV_Vect(Prism->Spline->Entry[number_of_splines].D, D);
1744
1745     if ((Prism->Spline_Type == QUADRATIC_SPLINE) ||
1746         (Prism->Spline_Type == CUBIC_SPLINE))
1747     {
1748       /* Get maximum coordinates in current segment. */
1749
1750       c[0] = 3.0 * A[X];
1751       c[1] = 2.0 * B[X];
1752       c[2] = C[X];
1753
1754       n = Solve_Polynomial(2, c, r, FALSE, 0.0);
1755
1756       while (n--)
1757       {
1758         if ((r[n] >= 0.0) && (r[n] <= 1.0))
1759         {
1760           x[n] = r[n] * (r[n] * (r[n] * A[X] + B[X]) + C[X]) + D[X];
1761         }
1762       }
1763
1764       c[0] = 3.0 * A[Y];
1765       c[1] = 2.0 * B[Y];
1766       c[2] = C[Y];
1767
1768       n = Solve_Polynomial(2, c, r, FALSE, 0.0);
1769
1770       while (n--)
1771       {
1772         if ((r[n] >= 0.0) && (r[n] <= 1.0))
1773         {
1774           y[n] = r[n] * (r[n] * (r[n] * A[Y] + B[Y]) + C[Y]) + D[Y];
1775         }
1776       }
1777     }
1778
1779     /* Set current segment's bounding rectangle. */
1780
1781     Prism->Spline->Entry[number_of_splines].x1 = min(min(x[0], x[1]), min(x[2], x[3]));
1782
1783     Prism->Spline->Entry[number_of_splines].x2 =
1784     Prism->Spline->Entry[number_of_splines].u2 = max(max(x[0], x[1]), max(x[2], x[3]));
1785
1786     Prism->Spline->Entry[number_of_splines].y1 =
1787     Prism->Spline->Entry[number_of_splines].v1 = min(min(y[0], y[1]), min(y[2], y[3]));
1788
1789     Prism->Spline->Entry[number_of_splines].y2 =
1790     Prism->Spline->Entry[number_of_splines].v2 = max(max(y[0], y[1]), max(y[2], y[3]));
1791
1792     /* Keep track of overall bounding rectangle. */
1793
1794     xmin = min(xmin, Prism->Spline->Entry[number_of_splines].x1);
1795     xmax = max(xmax, Prism->Spline->Entry[number_of_splines].x2);
1796
1797     ymin = min(ymin, Prism->Spline->Entry[number_of_splines].y1);
1798     ymax = max(ymax, Prism->Spline->Entry[number_of_splines].y2);
1799
1800     number_of_splines++;
1801
1802     /* Advance to next segment. */
1803
1804     switch (Prism->Spline_Type)
1805     {
1806       case LINEAR_SPLINE:
1807
1808         if ((fabs(P[i1][X] - First[X]) < EPSILON) &&
1809             (fabs(P[i1][Y] - First[Y]) < EPSILON))
1810         {
1811           i++;
1812
1813           if (i+1 < Prism->Number)
1814           {
1815             Assign_UV_Vect(First, P[i+1]);
1816           }
1817         }
1818
1819         break;
1820
1821       case QUADRATIC_SPLINE:
1822
1823         if ((fabs(P[i2][X] - First[X]) < EPSILON) &&
1824             (fabs(P[i2][Y] - First[Y]) < EPSILON))
1825         {
1826           i += 2;
1827
1828           if (i+2 < Prism->Number)
1829           {
1830             Assign_UV_Vect(First, P[i+2]);
1831           }
1832         }
1833
1834         break;
1835
1836       case CUBIC_SPLINE:
1837
1838         if ((fabs(P[i2][X] - First[X]) < EPSILON) &&
1839             (fabs(P[i2][Y] - First[Y]) < EPSILON))
1840         {
1841           i += 3;
1842
1843           if (i+2 < Prism->Number)
1844           {
1845             Assign_UV_Vect(First, P[i+2]);
1846           }
1847         }
1848
1849         break;
1850
1851       case BEZIER_SPLINE:
1852
1853         i += 3;
1854
1855         break;
1856     }
1857   }
1858
1859   Prism->Number = number_of_splines;
1860
1861   /* Set overall bounding rectangle. */
1862
1863   Prism->x1 =
1864   Prism->u1 = xmin;
1865   Prism->x2 =
1866   Prism->u2 = xmax;
1867
1868   Prism->y1 =
1869   Prism->v1 = ymin;
1870   Prism->y2 =
1871   Prism->v2 = ymax;
1872
1873   if (Prism->Sweep_Type == CONIC_SWEEP)
1874   {
1875     /* Recalculate bounding rectangles. */
1876
1877     for (i = 0; i < Prism->Number; i++)
1878     {
1879       x[0] = Prism->Spline->Entry[i].x1 * Prism->Height1;
1880       x[1] = Prism->Spline->Entry[i].x1 * Prism->Height2;
1881       x[2] = Prism->Spline->Entry[i].x2 * Prism->Height1;
1882       x[3] = Prism->Spline->Entry[i].x2 * Prism->Height2;
1883
1884       xmin = min(min(x[0], x[1]), min(x[2], x[3]));
1885       xmax = max(max(x[0], x[1]), max(x[2], x[3]));
1886
1887       Prism->Spline->Entry[i].x1 = xmin;
1888       Prism->Spline->Entry[i].x2 = xmax;
1889
1890       y[0] = Prism->Spline->Entry[i].y1 * Prism->Height1;
1891       y[1] = Prism->Spline->Entry[i].y1 * Prism->Height2;
1892       y[2] = Prism->Spline->Entry[i].y2 * Prism->Height1;
1893       y[3] = Prism->Spline->Entry[i].y2 * Prism->Height2;
1894
1895       ymin = min(min(y[0], y[1]), min(y[2], y[3]));
1896       ymax = max(max(y[0], y[1]), max(y[2], y[3]));
1897
1898       Prism->Spline->Entry[i].y1 = ymin;
1899       Prism->Spline->Entry[i].y2 = ymax;
1900     }
1901
1902     /* Recalculate overall bounding rectangle. */
1903
1904     x[0] = Prism->x1 * Prism->Height1;
1905     x[1] = Prism->x1 * Prism->Height2;
1906     x[2] = Prism->x2 * Prism->Height1;
1907     x[3] = Prism->x2 * Prism->Height2;
1908
1909     xmin = min(min(x[0], x[1]), min(x[2], x[3]));
1910     xmax = max(max(x[0], x[1]), max(x[2], x[3]));
1911
1912     Prism->x1 = xmin;
1913     Prism->x2 = xmax;
1914
1915     y[0] = Prism->y1 * Prism->Height1;
1916     y[1] = Prism->y1 * Prism->Height2;
1917     y[2] = Prism->y2 * Prism->Height1;
1918     y[3] = Prism->y2 * Prism->Height2;
1919
1920     ymin = min(min(y[0], y[1]), min(y[2], y[3]));
1921     ymax = max(max(y[0], y[1]), max(y[2], y[3]));
1922
1923     Prism->y1 = ymin;
1924     Prism->y2 = ymax;
1925   }
1926 }
1927
1928
1929