home *** CD-ROM | disk | FTP | other *** search
/ TOS Silver 2000 / TOS Silver 2000.iso / Anwendun / Pov / POVTOOLS / SOURCE / POLY.C < prev    next >
Encoding:
C/C++ Source or Header  |  1995-11-25  |  19.1 KB  |  759 lines
  1. /****************************************************************************
  2. *                poly.c
  3. *
  4. *  This module implements the code for general 3 variable polynomial shapes
  5. *
  6. *  This file was written by Alexander Enzmann.  He wrote the code for
  7. *  4th - 6th order shapes and generously provided us these enhancements.
  8. *
  9. *  from Persistence of Vision Raytracer
  10. *  Copyright 1993 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. If 
  17. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  18. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  19. *  Forum.  The latest version of POV-Ray may be found there as well.
  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. #include "frame.h"
  28. #include "vector.h"
  29. #include "povproto.h"
  30.  
  31. /* Basic form of a quartic equation
  32.    a00*x^4+a01*x^3*y+a02*x^3*z+a03*x^3+a04*x^2*y^2+
  33.    a05*x^2*y*z+a06*x^2*y+a07*x^2*z^2+a08*x^2*z+a09*x^2+
  34.    a10*x*y^3+a11*x*y^2*z+a12*x*y^2+a13*x*y*z^2+a14*x*y*z+
  35.    a15*x*y+a16*x*z^3+a17*x*z^2+a18*x*z+a19*x+a20*y^4+
  36.    a21*y^3*z+a22*y^3+a23*y^2*z^2+a24*y^2*z+a25*y^2+a26*y*z^3+
  37.    a27*y*z^2+a28*y*z+a29*y+a30*z^4+a31*z^3+a32*z^2+a33*z+a34
  38. */
  39.  
  40. #define POLYNOMIAL_TOLERANCE 1.0e-4
  41. #define COEFF_LIMIT 1.0e-20
  42. #define BINOMSIZE 40
  43.  
  44. /* The following table contains the binomial coefficients up to 15 */
  45. int binomials[15][15] =
  46.   {
  47.     {  
  48.     1,  0,  0,  0,   0,   0,   0,   0,   0,   0,   0,  0,  0,  0,  0
  49.   }
  50.   ,
  51.     {  
  52.     1,  1,  0,  0,   0,   0,   0,   0,   0,   0,   0,  0,  0,  0,  0
  53.   }
  54.   ,
  55.     {  
  56.     1,  2,  1,  0,   0,   0,   0,   0,   0,   0,   0,  0,  0,  0,  0
  57.   }
  58.   ,
  59.     {  
  60.     1,  3,  3,  1,   0,   0,   0,   0,   0,   0,   0,  0,  0,  0,  0
  61.   }
  62.   ,
  63.     {  
  64.     1,  4,  6,  4,   1,   0,   0,   0,   0,   0,   0,  0,  0,  0,  0
  65.   }
  66.   ,
  67.     {  
  68.     1,  5, 10, 10,   5,   1,   0,   0,   0,   0,   0,  0,  0,  0,  0
  69.   }
  70.   ,
  71.     {  
  72.     1,  6, 15, 20,  15,   6,   1,   0,   0,   0,   0,  0,  0,  0,  0
  73.   }
  74.   ,
  75.     {  
  76.     1,  7, 21, 35,  35,  21,   7,   1,   0,   0,   0,  0,  0,  0,  0
  77.   }
  78.   ,
  79.     {  
  80.     1,  8, 28, 56,  70,  56,  28,   8,   1,   0,   0,  0,  0,  0,  0
  81.   }
  82.   ,
  83.     {  
  84.     1,  9, 36, 84, 126, 126,  84,  36,   9,   1,   0,  0,  0,  0,  0
  85.   }
  86.   ,
  87.     {  
  88.     1, 10, 45,120, 210, 252, 210, 120,  45,  10,   1,  0,  0,  0,  0
  89.   }
  90.   ,
  91.     {  
  92.     1, 11, 55,165, 330, 462, 462, 330, 165,  55,  11,  1,  0,  0,  0
  93.   }
  94.   ,
  95.     {  
  96.     1, 12, 66,220, 495, 792, 924, 792, 495, 220,  66, 12,  1,  0,  0
  97.   }
  98.   ,
  99.     {  
  100.     1, 13, 78,286, 715,1287,1716,1716,1287, 715, 286, 78, 13,  1,  0
  101.   }
  102.   ,
  103.     {  
  104.     1, 14, 91,364,1001,2002,3003,3432,3003,2002,1001,364, 91, 14,  1
  105.   }
  106.   };
  107.  
  108. DBL eqn_v[3][MAX_ORDER+1], eqn_vt[3][MAX_ORDER+1];
  109.  
  110. METHODS Poly_Methods =
  111.   { 
  112.   All_Poly_Intersections,
  113.   Inside_Poly, Poly_Normal, Copy_Poly,
  114.   Translate_Poly, Rotate_Poly,
  115.   Scale_Poly, Transform_Poly, Invert_Poly, Destroy_Poly
  116. };
  117.  
  118. extern long Ray_Poly_Tests, Ray_Poly_Tests_Succeeded;
  119. extern unsigned int Options;
  120. extern int Shadow_Test_Flag;
  121.  
  122. /* unused
  123. static DBL evaluate_linear PARAMS((VECTOR *P, DBL *a));
  124. static DBL evaluate_quadratic PARAMS((VECTOR *P, DBL *a));
  125. */
  126.  
  127. static int intersect PARAMS((RAY *Ray, int Order, DBL *Coeffs, int Sturm_Flag,
  128. DBL *Depths));
  129. static void normal0 PARAMS((VECTOR *Result, int Order, DBL *Coeffs,
  130. VECTOR *IPoint));
  131. static void normal1 PARAMS((VECTOR *Result, int Order, DBL *Coeffs,
  132. VECTOR *IPoint));
  133. static DBL inside PARAMS((VECTOR *IPoint, int Order, DBL *Coeffs));
  134. static int intersect_linear PARAMS((RAY *ray, DBL *Coeffs, DBL *Depths));
  135. static int intersect_quadratic PARAMS((RAY *ray, DBL *Coeffs, DBL *Depths));
  136. static int factor_out PARAMS((int n, int i, int *c, int *s));
  137. static long binomial PARAMS((int n, int r));
  138. static void factor1 PARAMS((int n, int *c, int *s));
  139.  
  140. int All_Poly_Intersections(Object, Ray, Depth_Stack)
  141. OBJECT *Object;
  142. RAY *Ray;
  143. ISTACK *Depth_Stack;
  144.   {
  145.   POLY *Poly = (POLY *) Object;
  146.   DBL Depths[MAX_ORDER], len;
  147.   VECTOR IPoint, dv;
  148.   int cnt, i, j, Intersection_Found;
  149.   RAY New_Ray;
  150.  
  151.   /* Transform the ray into the polynomial's space */
  152.   if (Poly->Trans != NULL) 
  153.     {
  154.     MInvTransPoint(&New_Ray.Initial, &Ray->Initial, Poly->Trans);
  155.     MInvTransDirection(&New_Ray.Direction, &Ray->Direction, Poly->Trans);
  156.     }
  157.   else 
  158.     {
  159.     New_Ray.Initial   = Ray->Initial;
  160.     New_Ray.Direction = Ray->Direction;
  161.     }
  162.  
  163.   VDot(len, New_Ray.Direction, New_Ray.Direction);
  164.   if (len == 0.0)
  165.     return 0;
  166.   len = 1.0 / sqrt(len);
  167.   VScaleEq(New_Ray.Direction, len);
  168.  
  169.   Intersection_Found = FALSE;
  170.   Ray_Poly_Tests++;
  171.  
  172.   if (Poly->Order == 1)
  173.     cnt = intersect_linear(&New_Ray, Poly->Coeffs, Depths);
  174.   else if (Poly->Order == 2)
  175.     cnt = intersect_quadratic(&New_Ray, Poly->Coeffs, Depths);
  176.   else
  177.     cnt = intersect(&New_Ray, Poly->Order, Poly->Coeffs, Poly->Sturm_Flag,
  178.       Depths);
  179.   if (cnt > 0) Ray_Poly_Tests_Succeeded++;
  180.  
  181.   for (i=0;i<cnt;i++) 
  182.     {
  183.     if (Depths[i] < POLYNOMIAL_TOLERANCE) goto l0;
  184.     for (j=0;j<i;j++)
  185.       if (Depths[i] == Depths[j]) goto l0;
  186.     VScale(IPoint, New_Ray.Direction, Depths[i]);
  187.     VAddEq(IPoint, New_Ray.Initial);
  188.     /* Transform the point into world space */
  189.     if (Poly->Trans != NULL)
  190.       MTransPoint(&IPoint, &IPoint, Poly->Trans);
  191.  
  192.     VSub(dv, IPoint, Ray->Initial);
  193.     VLength(len, dv);
  194.     if (Point_In_Clip(&IPoint, Object->Clip))
  195.       {
  196.       push_entry(len,IPoint,Object,Depth_Stack);
  197.       Intersection_Found = TRUE;
  198.       }
  199. l0:;
  200.     }
  201.   return (Intersection_Found);
  202.   }
  203.  
  204. /* For speedup of low order polynomials, expand out the terms
  205.    involved in evaluating the poly. */
  206. /* unused
  207. static DBL
  208. evaluate_linear(P, a)
  209.    VECTOR *P;
  210.    DBL *a;
  211. {
  212.    return (a[0] * P->x) + (a[1] * P->y) + (a[2] * P->z) + a[3];
  213. }
  214.  
  215. static DBL
  216. evaluate_quadratic(P, a)
  217.    VECTOR *P;
  218.    DBL *a;
  219. {
  220.    DBL x, y, z;
  221.  
  222.    x = P->x; y = P->y; z = P->z;
  223.    return  a[0] * x * x + a[1] * x * y + a[2] * x * z +
  224.            a[3] * x     + a[4] * y * y + a[5] * y * z +
  225.            a[6] * y     + a[7] * z * z + a[8] * z     +
  226.            a[9];
  227. }
  228. */
  229. /* Remove all factors of i from n. */
  230. static int
  231. factor_out(n, i, c, s)
  232. int n, i, *c, *s;
  233.   {
  234.   while (!(n % i)) 
  235.     {
  236.     n /= i;
  237.     s[(*c)++] = i;
  238.     }
  239.   return n;
  240.   }
  241.  
  242. /* Find all prime factors of n. (Note that n must be less than 2^15 */
  243. static void
  244. factor1(n, c, s)
  245. int n, *c, *s;
  246.   {
  247.   int i,k;
  248.   /* First factor out any 2s */
  249.   n = factor_out(n, 2, c, s);
  250.   /* Now any odd factors */
  251.   k = (int)sqrt(n) + 1;
  252.   for (i=3;n>1 && i<=k;i+=2)
  253.     if (!(n%i)) 
  254.     {
  255.     n = factor_out(n, i, c, s);
  256.     k = (int)sqrt(n)+1;
  257.     }
  258.   if (n>1)
  259.     s[(*c)++] = n;
  260.   }
  261.  
  262. /* Calculate the binomial coefficent of n,r. */
  263. static
  264. long
  265. binomial(n, r)
  266. int n, r;
  267.   {
  268.   int h,i,j,k,l;
  269.   unsigned long result;
  270.   static int stack1[BINOMSIZE], stack2[BINOMSIZE];
  271.   if (n<0 || r<0 || r>n)
  272.     result = 0L;
  273.   else if (r==n)
  274.     result = 1L;
  275.   else if (r < 15 && n < 15)
  276.     result = (long)binomials[n][r];
  277.   else 
  278.     {
  279.     j = 0;
  280.     for (i=r+1;i<=n;i++)
  281.       stack1[j++] = i;
  282.     for (i=2;i<=(n-r);i++) 
  283.       {
  284.       h = 0;
  285.       factor1(i, &h, stack2);
  286.       for (k=0;k<h;k++) 
  287.         {
  288.         for (l=0;l<j;l++)
  289.           if (!(stack1[l] % stack2[k])) 
  290.           {
  291.           stack1[l] /= stack2[k];
  292.           goto l1;
  293.           }
  294.         /* Error if we get here */
  295.         if (Options & DEBUGGING) 
  296.           {
  297.           printf("Failed to factor\n");
  298.           fflush(stdout);
  299.           }
  300. l1:;
  301.         }
  302.       }
  303.     result=1;
  304.     for (i=0;i<j;i++)
  305.       result *= stack1[i];
  306.     }
  307.   return result;
  308.   }
  309.  
  310. static DBL
  311. inside(IPoint, Order, Coeffs)
  312. VECTOR *IPoint;
  313. int Order;
  314. DBL *Coeffs;
  315.   {
  316.   DBL x[MAX_ORDER+1], y[MAX_ORDER+1];
  317.   DBL z[MAX_ORDER+1], c, Result;
  318.   int i, j, k, term;
  319.  
  320.   x[0] = 1.0;       y[0] = 1.0;       z[0] = 1.0;
  321.   x[1] = IPoint->x; y[1] = IPoint->y; z[1] = IPoint->z;
  322.   for (i=2;i<=Order;i++) 
  323.     {
  324.     x[i] = x[1] * x[i-1];
  325.     y[i] = y[1] * y[i-1];
  326.     z[i] = z[1] * z[i-1];
  327.     }
  328.   Result = 0.0;
  329.   term = 0;
  330.   for (i=Order;i>=0;i--)
  331.     for (j=Order-i;j>=0;j--)
  332.     for (k=Order-(i+j);k>=0;k--) 
  333.     {
  334.     if ((c = Coeffs[term]) != 0.0)
  335.       Result += c * x[i] * y[j] * z[k];
  336.     term++;
  337.     }
  338.   return Result;
  339.   }
  340.  
  341. /* Intersection of a ray and an arbitrary polynomial function */
  342. static int
  343. intersect(ray, Order, Coeffs, Sturm_Flag, Depths)
  344. RAY *ray;
  345. int Order, Sturm_Flag;
  346. DBL *Coeffs, *Depths;
  347.   {
  348.   DBL eqn[MAX_ORDER+1];
  349.   DBL t[3][MAX_ORDER+1];
  350.   VECTOR P, D;
  351.   DBL val;
  352.   int h, i, j, k, i1, j1, k1, term;
  353.   int offset;
  354.  
  355.   /* First we calculate the values of the individual powers
  356.       of x, y, and z as they are represented by the ray */
  357.   P = ray->Initial;
  358.   D = ray->Direction;
  359.   for (i=0;i<3;i++) 
  360.     {
  361.     eqn_v[i][0]  = 1.0;
  362.     eqn_vt[i][0] = 1.0;
  363.     }
  364.   eqn_v[0][1] = P.x;
  365.   eqn_v[1][1] = P.y;
  366.   eqn_v[2][1] = P.z;
  367.   eqn_vt[0][1] = D.x;
  368.   eqn_vt[1][1] = D.y;
  369.   eqn_vt[2][1] = D.z;
  370.  
  371.   for (i=2;i<=Order;i++)
  372.     for (j=0;j<3;j++) 
  373.     {
  374.     eqn_v[j][i]  = eqn_v[j][1] * eqn_v[j][i-1];
  375.     eqn_vt[j][i] = eqn_vt[j][1] * eqn_vt[j][i-1];
  376.     }
  377.  
  378.   for (i=0;i<=Order;i++)
  379.     eqn[i] = 0.0;
  380.  
  381.   /* Now walk through the terms of the polynomial.  As we go
  382.       we substitute the ray equation for each of the variables. */
  383.   term = 0;
  384.   for (i=Order;i>=0;i--) 
  385.     {
  386.     for (h=0;h<=i;h++)
  387.       t[0][h] = binomial(i, h) *
  388.       eqn_vt[0][i-h] *
  389.       eqn_v[0][h];
  390.     for (j=Order-i;j>=0;j--) 
  391.       {
  392.       for (h=0;h<=j;h++)
  393.         t[1][h] = binomial(j, h) *
  394.         eqn_vt[1][j-h] *
  395.         eqn_v[1][h];
  396.       for (k=Order-(i+j);k>=0;k--) 
  397.         {
  398.         if (Coeffs[term] != 0) 
  399.           {
  400.           for (h=0;h<=k;h++)
  401.             t[2][h] = binomial(k, h) *
  402.             eqn_vt[2][k-h] *
  403.             eqn_v[2][h];
  404.  
  405.           /* Multiply together the three polynomials */
  406.           offset = Order - (i + j + k);
  407.           for (i1=0;i1<=i;i1++)
  408.             for (j1=0;j1<=j;j1++)
  409.             for (k1=0;k1<=k;k1++) 
  410.             {
  411.             h = offset + i1 + j1 + k1;
  412.             val = Coeffs[term];
  413.             val *= t[0][i1];
  414.             val *= t[1][j1];
  415.             val *= t[2][k1];
  416.             eqn[h] += val;
  417.             }
  418.           }
  419.         term++;
  420.         }
  421.       }
  422.     }
  423.   for (i=0,j=Order;i<=Order;i++)
  424.     if (eqn[i] != 0.0)
  425.       break;
  426.     else
  427.       j--;
  428.  
  429.   if (j > 4 || Sturm_Flag)
  430.     return polysolve(j, &eqn[i], Depths);
  431.   else if (j == 4)
  432.     return solve_quartic(&eqn[i], Depths);
  433.   else if (j==3)
  434.     return solve_cubic(&eqn[i], Depths);
  435.   else if (j==2)
  436.     return solve_quadratic(&eqn[i], Depths);
  437.   else
  438.     return 0;
  439.   }
  440.  
  441. /* Intersection of a ray and a quadratic */
  442. static int
  443. intersect_linear(ray, Coeffs, Depths)
  444. RAY *ray;
  445. DBL *Coeffs, *Depths;
  446.   {
  447.   DBL t0, t1, *a = Coeffs;
  448.  
  449.   t0 = a[0] * ray->Initial.x + a[1] * ray->Initial.y + a[2] * ray->Initial.z;
  450.   t1 = a[0] * ray->Direction.x + a[1] * ray->Direction.y +
  451.   a[2] * ray->Direction.z;
  452.   if (fabs(t1) < EPSILON)
  453.     return 0;
  454.   Depths[0] = -(a[3] + t0) / t1;
  455.   return 1;
  456.   }
  457.  
  458. /* Intersection of a ray and a quadratic */
  459. static int
  460. intersect_quadratic(ray, Coeffs, Depths)
  461. RAY *ray;
  462. DBL *Coeffs, *Depths;
  463.   {
  464.   DBL x, y, z, x2, y2, z2;
  465.   DBL xx, yy, zz, xx2, yy2, zz2;
  466.   DBL *a, ac, bc, cc, d, t;
  467.  
  468.   x  = ray->Initial.x;
  469.   y  = ray->Initial.y;
  470.   z  = ray->Initial.z;
  471.   xx = ray->Direction.x;
  472.   yy = ray->Direction.y;
  473.   zz = ray->Direction.z;
  474.   x2 = x * x;  y2 = y * y;  z2 = z * z;
  475.   xx2 = xx * xx;  yy2 = yy * yy;  zz2 = zz * zz;
  476.   a = Coeffs;
  477.  
  478.   /*
  479.       Determine the coefficients of t^n, where the line is represented
  480.       as (x,y,z) + (xx,yy,zz)*t.
  481.    */
  482.   ac = (a[0]*xx2 + a[1]*xx*yy + a[2]*xx*zz + a[4]*yy2 + a[5]*yy*zz +
  483.     a[7]*zz2);
  484.   bc = (2*a[0]*x*xx + a[1]*(x*yy + xx*y) + a[2]*(x*zz + xx*z) +
  485.     a[3]*xx + 2*a[4]*y*yy + a[5]*(y*zz + yy*z) + a[6]*yy +
  486.     2*a[7]*z*zz + a[8]*zz);
  487.   cc = a[0]*x2 + a[1]*x*y + a[2]*x*z + a[3]*x + a[4]*y2 +
  488.   a[5]*y*z + a[6]*y + a[7]*z2 + a[8]*z + a[9];
  489.  
  490.   if (fabs(ac) < COEFF_LIMIT) 
  491.     {
  492.     if (fabs(bc) < COEFF_LIMIT)
  493.       return 0;
  494.     Depths[0] = -cc / bc;
  495.     return 1;
  496.     }
  497.  
  498.   /* Solve the quadratic formula & return results that are
  499.       within the correct interval. */
  500.   d = bc * bc - 4.0 * ac * cc;
  501.   if (d < 0.0) return 0;
  502.   d = sqrt(d);
  503.   bc = -bc;
  504.   t = 2.0 * ac;
  505.   Depths[0] = (bc + d) / t;
  506.   Depths[1] = (bc - d) / t;
  507.   return 2;
  508.   }
  509.  
  510. /* Normal to a polynomial (used for polynomials with order > 4) */
  511. static void normal0(Result, Order, Coeffs, IPoint)
  512. VECTOR *Result;
  513. int Order;
  514. DBL *Coeffs;
  515. VECTOR *IPoint;
  516.   {
  517.   int i, j, k, term;
  518.   DBL val, *a, x[MAX_ORDER+1], y[MAX_ORDER+1], z[MAX_ORDER+1];
  519.  
  520.   x[0] = 1.0; y[0] = 1.0; z[0] = 1.0;
  521.   x[1] = IPoint->x;
  522.   y[1] = IPoint->y;
  523.   z[1] = IPoint->z;
  524.   for (i=2;i<=Order;i++) 
  525.     {
  526.     x[i] = IPoint->x * x[i-1];
  527.     y[i] = IPoint->y * y[i-1];
  528.     z[i] = IPoint->z * z[i-1];
  529.     }
  530.   a = Coeffs;
  531.   term = 0;
  532.   Make_Vector(Result, 0, 0, 0);
  533.   for (i=Order;i>=0;i--)
  534.     for (j=Order-i;j>=0;j--)
  535.     for (k=Order-(i+j);k>=0;k--) 
  536.     {
  537.     if (i >= 1) 
  538.       {
  539.       val = x[i-1] * y[j] * z[k];
  540.       Result->x += i * a[term] * val;
  541.       }
  542.     if (j >= 1) 
  543.       {
  544.       val = x[i] * y[j-1] * z[k];
  545.       Result->y += j * a[term] * val;
  546.       }
  547.     if (k >= 1) 
  548.       {
  549.       val = x[i] * y[j] * z[k-1];
  550.       Result->z += k * a[term] * val;
  551.       }
  552.     term++;
  553.     }
  554.   }
  555.  
  556. /* Normal to a polynomial (for polynomials of order <= 4) */
  557. static void
  558. normal1(Result, Order, Coeffs, IPoint)
  559. VECTOR *Result;
  560. int Order;
  561. DBL *Coeffs;
  562. VECTOR *IPoint;
  563.   {
  564.   DBL *a, x, y, z, x2, y2, z2, x3, y3, z3;
  565.  
  566.   a = Coeffs;
  567.   x = IPoint->x;
  568.   y = IPoint->y;
  569.   z = IPoint->z;
  570.  
  571.   if (Order == 1)
  572.     /* Linear partial derivatives */
  573.     Make_Vector(Result, a[0], a[1], a[2])
  574. else if (Order == 2) 
  575.   {
  576.   /* Quadratic partial derivatives */
  577.   Result->x = 2*a[0]*x+a[1]*y+a[2]*z+a[3];
  578.   Result->y = a[1]*x+2*a[4]*y+a[5]*z+a[6];
  579.   Result->z = a[2]*x+a[5]*y+2*a[7]*z+a[8];
  580.   }
  581. else if (Order == 3) 
  582.   {
  583.   x2 = x * x;  y2 = y * y;  z2 = z * z;
  584.   /* Cubic partial derivatives */
  585.   Result->x = 3*a[0]*x2 + 2*x*(a[1]*y + a[2]*z + a[3]) + a[4]*y2 +
  586.   y*(a[5]*z + a[6]) + a[7]*z2 + a[8]*z + a[9];
  587.   Result->y = a[1]*x2 + x*(2*a[4]*y + a[5]*z + a[6]) + 3*a[10]*y2 +
  588.   2*y*(a[11]*z + a[12]) + a[13]*z2 + a[14]*z + a[15];
  589.   Result->z = a[2]*x2 + x*(a[5]*y + 2*a[7]*z + a[8]) + a[11]*y2 +
  590.   y*(2*a[13]*z + a[14]) + 3*a[16]*z2 + 2*a[17]*z + a[18];
  591.   }
  592. else 
  593.   {
  594.   /* Quartic partial derivatives */
  595.     x2 = x * x;  y2 = y * y;  z2 = z * z;
  596.   x3 = x * x2; y3 = y * y2; z3 = z * z2;
  597.   Result->x = 4*a[ 0]*x3+3*x2*(a[ 1]*y+a[ 2]*z+a[ 3])+
  598.   2*x*(a[ 4]*y2+y*(a[ 5]*z+a[ 6])+a[ 7]*z2+a[ 8]*z+a[ 9])+
  599.   a[10]*y3+y2*(a[11]*z+a[12])+y*(a[13]*z2+a[14]*z+a[15])+
  600.   a[16]*z3+a[17]*z2+a[18]*z+a[19];
  601.   Result->y = a[ 1]*x3+x2*(2*a[ 4]*y+a[ 5]*z+a[ 6])+
  602.   x*(3*a[10]*y2+2*y*(a[11]*z+a[12])+a[13]*z2+a[14]*z+a[15])+
  603.   4*a[20]*y3+3*y2*(a[21]*z+a[22])+2*y*(a[23]*z2+a[24]*z+a[25])+
  604.   a[26]*z3+a[27]*z2+a[28]*z+a[29];
  605.   Result->z = a[ 2]*x3+x2*(a[ 5]*y+2*a[ 7]*z+a[ 8])+
  606.   x*(a[11]*y2+y*(2*a[13]*z+a[14])+3*a[16]*z2+2*a[17]*z+a[18])+
  607.   a[21]*y3+y2*(2*a[23]*z+a[24])+y*(3*a[26]*z2+2*a[27]*z+a[28])+
  608.   4*a[30]*z3+3*a[31]*z2+2*a[32]*z+a[33];
  609.   }
  610. }
  611.  
  612. int Inside_Poly (IPoint, Object)
  613. VECTOR *IPoint;
  614. OBJECT *Object;
  615. {
  616. VECTOR New_Point;
  617. DBL Result;
  618.  
  619. /* Transform the point into polynomial's space */
  620. if (((POLY *)Object)->Trans != NULL)
  621. MInvTransPoint(&New_Point, IPoint, ((POLY *)Object)->Trans);
  622. else
  623.   New_Point = *IPoint;
  624.  
  625. Result = inside(&New_Point, ((POLY *)Object)->Order, ((POLY *)Object)->Coeffs);
  626. if (Result < Small_Tolerance)
  627. return ((int)(1-((POLY *)Object)->Inverted));
  628. else
  629.   return ((int)((POLY *)Object)->Inverted);
  630. }
  631.  
  632. /* Normal to a polynomial */
  633. void Poly_Normal(Result, Object, IPoint)
  634. VECTOR *Result, *IPoint;
  635. OBJECT *Object;
  636. {
  637. DBL val;
  638. VECTOR New_Point;
  639. POLY *Shape = (POLY *)Object;
  640.  
  641. /* Transform the point into the polynomials space */
  642. if (Shape->Trans != NULL)
  643. MInvTransPoint(&New_Point, IPoint, Shape->Trans);
  644. else
  645.   New_Point = *IPoint;
  646.  
  647. if (((POLY *)Object)->Order > 4)
  648. normal0(Result, Shape->Order, Shape->Coeffs, &New_Point);
  649. else
  650.   normal1(Result, Shape->Order, Shape->Coeffs, &New_Point);
  651.  
  652. /* Transform back to world space */
  653. if (Shape->Trans != NULL)
  654. MTransNormal(Result, Result, Shape->Trans);
  655.  
  656. /* Normalize (accounting for the possibility of a 0 length normal) */
  657. VDot(val, *Result, *Result);
  658. if (val > 0.0) 
  659.   {
  660.   val = 1.0 / sqrt(val);
  661.   VScaleEq(*Result, val);
  662.   }
  663. else
  664.   Make_Vector(Result, 1, 0, 0)
  665.     }
  666.  
  667.   /* Make a copy of a polynomial object */
  668.   void *Copy_Poly(Object)
  669.     OBJECT *Object;
  670. {
  671. POLY *New;
  672. int i;
  673.  
  674. New = Create_Poly (((POLY *)Object)->Order);
  675. New->Sturm_Flag = ((POLY *)Object)->Sturm_Flag;
  676. New->Inverted = ((POLY *)Object)->Inverted;
  677.  
  678. New->Trans = Copy_Transform(((POLY *)Object)->Trans);
  679.  
  680. for (i=0;i<term_counts(New->Order);i++)
  681. New->Coeffs[i] = ((POLY *)Object)->Coeffs[i];
  682. return (New);
  683. }
  684.  
  685. void Translate_Poly (Object, Vector)
  686. OBJECT *Object;
  687. VECTOR *Vector;
  688. {
  689. TRANSFORM Trans;
  690. Compute_Translation_Transform(&Trans, Vector);
  691. Transform_Poly(Object, &Trans);
  692. }
  693.  
  694. void Rotate_Poly (Object, Vector)
  695. OBJECT *Object;
  696. VECTOR *Vector;
  697. {
  698. TRANSFORM Trans;
  699. Compute_Rotation_Transform(&Trans, Vector);
  700. Transform_Poly(Object, &Trans);
  701. }
  702.  
  703. void Scale_Poly (Object, Vector)
  704. OBJECT *Object;
  705. VECTOR *Vector;
  706. {
  707. TRANSFORM Trans;
  708. Compute_Scaling_Transform(&Trans, Vector);
  709. Transform_Poly(Object, &Trans);
  710. }
  711.  
  712. void Transform_Poly(Object,Trans)
  713. OBJECT *Object;
  714. TRANSFORM *Trans;
  715. {
  716. if (((POLY *)Object)->Trans == NULL)
  717. ((POLY *)Object)->Trans = Create_Transform();
  718. recompute_bbox(&Object->Bounds, Trans);
  719. Compose_Transforms(((POLY *)Object)->Trans, Trans);
  720. }
  721.  
  722. void Invert_Poly(Object)
  723. OBJECT *Object;
  724. {
  725. ((POLY *) Object)->Inverted = 1 - ((POLY *)Object)->Inverted;
  726. }
  727.  
  728. POLY *Create_Poly(Order)
  729. int Order;
  730. {
  731. POLY *New;
  732. int i;
  733.  
  734. if ((New = (POLY *) malloc (sizeof (POLY))) == NULL)
  735. MAError ("poly");
  736.  
  737. INIT_OBJECT_FIELDS(New,POLY_OBJECT, &Poly_Methods);
  738. New->Order = Order;
  739. New->Sturm_Flag = FALSE;
  740. New->Trans = NULL;
  741. New->Inverted = FALSE;
  742. New->Coeffs = (DBL *)malloc(term_counts(Order) * sizeof(DBL));
  743. if (New->Coeffs == NULL)
  744. MAError("coefficients for POLY");
  745. for (i=0;i<term_counts(Order);i++)
  746. New->Coeffs[i] = 0.0;
  747.  
  748. return (New);
  749. }
  750.  
  751. void Destroy_Poly(Object)
  752. OBJECT *Object;
  753. {
  754. Destroy_Transform (((POLY *)Object)->Trans);
  755. free (((POLY *)Object)->Coeffs);
  756. free (Object);
  757. }
  758.  
  759.