home *** CD-ROM | disk | FTP | other *** search
/ Otherware / Otherware_1_SB_Development.iso / mac / developm / source / povsrc.sit / SOURCE / POLY.C < prev    next >
Encoding:
C/C++ Source or Header  |  1992-07-03  |  25.0 KB  |  819 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 1992 Persistence of Vision Team
  11. *---------------------------------------------------------------------------
  12. *  Copying, distribution and legal info is in the file povlegal.doc which
  13. *  should be distributed with this file. If povlegal.doc is not available
  14. *  or for more info please contact:
  15. *
  16. *       Drew Wells [POV-Team Leader] 
  17. *       CIS: 73767,1244  Internet: 73767.1244@compuserve.com
  18. *       Phone: (213) 254-4041
  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. #include "frame.h"
  27. #include "vector.h"
  28. #include "povproto.h"
  29.  
  30. /* Basic form of a quartic equation
  31.    a00*x^4+a01*x^3*y+a02*x^3*z+a03*x^3+a04*x^2*y^2+
  32.    a05*x^2*y*z+a06*x^2*y+a07*x^2*z^2+a08*x^2*z+a09*x^2+
  33.    a10*x*y^3+a11*x*y^2*z+a12*x*y^2+a13*x*y*z^2+a14*x*y*z+
  34.    a15*x*y+a16*x*z^3+a17*x*z^2+a18*x*z+a19*x+a20*y^4+
  35.    a21*y^3*z+a22*y^3+a23*y^2*z^2+a24*y^2*z+a25*y^2+a26*y*z^3+
  36.    a27*y*z^2+a28*y*z+a29*y+a30*z^4+a31*z^3+a32*z^2+a33*z+a34
  37. */
  38.  
  39. #define COEFF_LIMIT 1.0e-20
  40.  
  41. #define LEFT_MARGIN 0
  42. #define RIGHT_MARGIN 72
  43.  
  44. int binomial[11][12] = {
  45.    { 0, 1,  0,  0,   0,   0,   0,   0,   0, 0,  0,  0},
  46.    { 0, 1,  1,  0,   0,   0,   0,   0,   0, 0,  0,  0},
  47.    { 0, 1,  2,  1,   0,   0,   0,   0,   0, 0,  0,  0},
  48.    { 0, 1,  3,  3,   1,   0,   0,   0,   0, 0,  0,  0},
  49.    { 0, 1,  4,  6,   4,   1,   0,   0,   0, 0,  0,  0},
  50.    { 0, 1,  5, 10,  10,   5,   1,   0,   0, 0,  0,  0},
  51.    { 0, 1,  6, 15,  20,  15,   6,   1,   0, 0,  0,  0},
  52.    { 0, 1,  7, 21,  35,  35,  21,   7,   1, 0,  0,  0},
  53.    { 0, 1,  8, 28,  56,  70,  56,  28,   8, 1,  0,  0},
  54.    { 0, 1,  9, 36,  84, 126, 126,  84,  36, 9,  1,  0},
  55.    { 0, 1, 10, 45, 120, 210, 252, 210, 120, 45, 10, 1}};
  56.  
  57. int factorials[MAX_ORDER+1] = {1, 1, 2, 6, 24, 120, 720, 5040};
  58. int term_counts[MAX_ORDER+1] = {1, 4, 10, 20, 35, 56, 84, 120};
  59.  
  60. METHODS Poly_Methods =
  61. { Object_Intersect, All_Poly_Intersections,
  62.    Inside_Poly, Poly_Normal, Copy_Poly,
  63.    Translate_Poly, Rotate_Poly,
  64.    Scale_Poly, Invert_Poly};
  65.  
  66. extern long Ray_Poly_Tests, Ray_Poly_Tests_Succeeded;
  67.  
  68. static void transform PARAMS((int order,DBL *Coeffs,MATRIX *q));
  69. static int roll PARAMS((int order,int x,int y,int z));
  70. static void unroll PARAMS((int order,int index,int *x,int *y,int *z,int *w));
  71. /*static void show_poly PARAMS((char *label, int Order, DBL *Coeffs));*/
  72. static int intersect PARAMS((RAY *Ray, int Order, DBL *Coeffs, DBL *Depths));
  73. static int intersect_quartic PARAMS((RAY *Ray, POLY *Shape, DBL *Depths));
  74. static void quartic_normal PARAMS((VECTOR *Result, OBJECT *Object,
  75. VECTOR *Intersection_Point));
  76. static DBL inside PARAMS((VECTOR *point, int Order, DBL *Coeffs));
  77. static void normalp PARAMS((VECTOR *Result, int Order, DBL *Coeffs,
  78. VECTOR *Intersection_Point));
  79. static DBL do_partial_term PARAMS((MATRIX *q, int row, int pwr, int i,
  80. int j, int k, int l));
  81.  
  82. int All_Poly_Intersections(Object, Ray, Depth_Queue)
  83. OBJECT *Object;
  84. RAY *Ray;
  85. PRIOQ *Depth_Queue;
  86. {
  87.    POLY *Shape = (POLY *) Object;
  88.    DBL Depths[MAX_ORDER], len;
  89.    VECTOR Intersection_Point, dv;
  90.    INTERSECTION Local_Element;
  91.    int cnt, i, j, Intersection_Found;
  92.    RAY New_Ray;
  93.  
  94.    /* Transform the ray into the polynomial's space */
  95.    if (Shape->Transform != NULL) {
  96.       MInverseTransformVector(&New_Ray.Initial, &Ray->Initial,
  97.          Shape->Transform);
  98.       MInvTransVector(&New_Ray.Direction, &Ray->Direction, Shape->Transform);
  99.    }
  100.    else {
  101.       New_Ray.Initial.x = Ray->Initial.x;
  102.       New_Ray.Initial.y = Ray->Initial.y;
  103.       New_Ray.Initial.z = Ray->Initial.z;
  104.       New_Ray.Direction.x = Ray->Direction.x;
  105.       New_Ray.Direction.y = Ray->Direction.y;
  106.       New_Ray.Direction.z = Ray->Direction.z;
  107.    }
  108.  
  109.       len = sqrt(New_Ray.Direction.x * New_Ray.Direction.x +
  110.          New_Ray.Direction.y * New_Ray.Direction.y +
  111.          New_Ray.Direction.z * New_Ray.Direction.z);
  112.    if (len == 0.0)
  113.       return 0;
  114.    else {
  115.       New_Ray.Direction.x /= len;
  116.       New_Ray.Direction.y /= len;
  117.       New_Ray.Direction.z /= len;
  118.    }
  119.  
  120.       Intersection_Found = FALSE;
  121.    Ray_Poly_Tests++;
  122.    if (Shape->Order == 4)
  123.       cnt = intersect_quartic(&New_Ray, Shape, Depths);
  124.    else
  125.       cnt = intersect(&New_Ray, Shape->Order, Shape->Coeffs, &Depths[0]);
  126.  
  127.    if (cnt > 0) Ray_Poly_Tests_Succeeded++;
  128.    for (i=0;i<cnt;i++) {
  129.       if (Depths[i] < 0) goto l0;
  130.       for (j=0;j<i;j++)
  131.          if (Depths[i] == Depths[j]) goto l0;
  132.       VScale(Intersection_Point, New_Ray.Direction, Depths[j]);
  133.       VAdd(Intersection_Point, Intersection_Point, New_Ray.Initial);
  134.       /* Transform the point into world space */
  135.       if (Shape->Transform != NULL)
  136.          MTransformVector(&Intersection_Point, &Intersection_Point,
  137.             Shape->Transform);
  138.  
  139.       VSub(dv, Intersection_Point, Ray->Initial);
  140.       VLength(len, dv);
  141.       Local_Element.Depth = len;
  142.       Local_Element.Object = Shape->Parent_Object;
  143.       Local_Element.Point = Intersection_Point;
  144.       Local_Element.Shape = (SHAPE *)Shape;
  145.       pq_add (Depth_Queue, &Local_Element);
  146.       Intersection_Found = TRUE;
  147. l0:;
  148.    }
  149.    return (Intersection_Found);
  150. }
  151.  
  152. /* Given the powers return the index into the polynomial */
  153. static int roll(order, x, y, z)
  154. int order, x, y, z;
  155. {
  156.    int xstart, ystart, zstart;
  157.    xstart = binomial[order-x+2][order-x];
  158.    order = order - x;
  159.    ystart = binomial[order-y+1][order-y];
  160.    order = order - y;
  161.    zstart = binomial[order-z][order-z];
  162.    return xstart+ystart+zstart;
  163. }
  164.  
  165. /* Given the index into the polynomial, return the powers. */
  166. static void unroll(order, index, x, y, z, w)
  167. int order, index, *x, *y, *z, *w;
  168. {
  169.    int i, torder;
  170.    if (index==0) { *x = order; *y = 0; *z = 0; *w = 0; return; }
  171.    else if (order==1)
  172.       if (index == 1) { *x = 0; *y = 1; *z = 0; *w = 0; return; }
  173.       else { *x = 0; *y = 0; *z = 3 - index; *w = order-(*x+*y+*z); return; }
  174.    else for (i=0;binomial[3+i][i+1]<=index;i++);
  175.    torder = order;
  176.    *x = torder - i; index -= binomial[2+i][i]; torder = i;
  177.    if (index==0) { *y = torder; *z = 0; *w = 0; return; }
  178.    else if (torder==1) { *y = 0; *z = 2-index; *w = order-(*x+*y+*z); return; }
  179.    else for (i=0;binomial[2+i][i+1]<=index;i++);
  180.    *y = torder - i; index -= binomial[1+i][i]; torder = i;
  181.    *z = torder - index;
  182.    *w = order - (*x + *y + *z);
  183. }
  184.  
  185.    /* REMOVED UNUSED ROUTINE 3/27/92
  186. static void show_poly(label, Order, Coeffs)
  187.    char *label;
  188.    int Order;
  189.    DBL *Coeffs;
  190. {
  191.    int i,j,cnt,col;
  192.    int wp, xp, yp, zp;
  193.    char term[40];
  194.    j = 0;
  195.    printf("%s", label);
  196.    col = strlen(label);
  197.  
  198.  
  199.    for (i=0;i<term_counts[Order];i++)
  200.       if (Coeffs[i] != 0.0) {
  201.      if (j) cnt = STRLN(sprintf(&term[0], " + %.5lg", Coeffs[i]));
  202.      else cnt = STRLN(sprintf(&term[0], "%.5lg", Coeffs[i]));
  203.      unroll(Order, i, &xp, &yp, &zp, &wp);
  204.      if (xp == 1) cnt += STRLN(sprintf(&term[cnt], "*x"));
  205.      else if (xp > 1) cnt += STRLN(sprintf(&term[cnt], "*x^%d", xp));
  206.      if (yp == 1) cnt += STRLN(sprintf(&term[cnt], "*y"));
  207.      else if (yp > 1) cnt += STRLN(sprintf(&term[cnt], "*y^%d", yp));
  208.      if (zp == 1) cnt += STRLN(sprintf(&term[cnt], "*z"));
  209.      else if (zp > 1) cnt += STRLN(sprintf(&term[cnt], "*z^%d", zp));
  210.      j = 1;
  211.      if (col+cnt>RIGHT_MARGIN) { printf("\n"); col=LEFT_MARGIN; }
  212.      printf("%s", &term[0]);
  213.      col += cnt;
  214.      }
  215. }
  216. */
  217.  
  218. /* Intersection of a ray and an arbitrary polynomial function */
  219. static int intersect(Ray, Order, Coeffs, Depths)
  220. RAY *Ray;
  221. int Order;
  222. DBL *Coeffs, *Depths;
  223. {
  224.    MATRIX q;
  225.    DBL *a, t[MAX_ORDER+1];
  226.    int i, j;
  227.    /* Determine the coefficients of t^n, where the line is represented
  228.       as (x,y,z) + (xx,yy,zz)*t.  */
  229.    if ((a = (DBL *)malloc(term_counts[Order] * sizeof(DBL))) == NULL) {
  230.       printf("Cannot allocate memory for coefficients in poly intersect()\n");
  231.       exit(1);
  232.    }
  233.    for (i=0;i<term_counts[Order];i++)
  234.       a[i] = Coeffs[i];
  235.    MZero((MATRIX *)&q[0][0]);
  236.    q[0][0] = Ray->Direction.x; q[3][0] = Ray->Initial.x;
  237.    q[0][1] = Ray->Direction.y; q[3][1] = Ray->Initial.y;
  238.    q[0][2] = Ray->Direction.z; q[3][2] = Ray->Initial.z;
  239.    transform(Order, a, (MATRIX *)&q[0][0]);
  240.    /* The equation is now in terms of one variable.  Use numerical
  241.       techniques to solve the polynomial that represents the intersections. */
  242.    for (i=0;i<=Order;i++) {
  243.       t[i] = a[binomial[3+i][4]-1];
  244.       if (t[i] > -COEFF_LIMIT && t[i] < COEFF_LIMIT)
  245.          t[i] = 0.0;
  246.    }
  247.    free(a);
  248.    for (i=0,j=Order;i<=Order;i++)
  249.       if (t[i] != 0.0)
  250.          break;
  251.       else
  252.          j -= 1;
  253.    if (j > 2)
  254.       return polysolve(j, &t[i], Depths);
  255.    else if (j > 0)
  256.       return solve_quadratic(&t[i], Depths);
  257.    else
  258.       return 0;
  259. }
  260.  
  261. static DBL inside(point, Order, Coeffs)
  262. VECTOR *point;
  263. int Order;
  264. DBL *Coeffs;
  265. {
  266.    DBL x[MAX_ORDER+1], y[MAX_ORDER+1], z[MAX_ORDER+1], Result;
  267.    int i, k0, k1, k2, k3;
  268.    x[0] = 1.0; y[0] = 1.0; z[0] = 1.0;
  269.    x[1] = point->x; y[1] = point->y; z[1] = point->z;
  270.    for (i=2;i<=MAX_ORDER;i++) {
  271.       x[i] = x[1] * x[i-1]; y[i] = y[1] * y[i-1]; z[i] = z[1] * z[i-1];
  272.    }
  273.    Result = 0.0;
  274.    for (i=0;i<term_counts[Order];i++) {
  275.       unroll(Order, i, &k0, &k1, &k2, &k3);
  276.       Result += Coeffs[i] * x[k0] * y[k1] * z[k2];
  277.    }
  278.  
  279.    /* The Epsilon fudge factor is so that points really near the
  280.       surface are considered inside the surface */
  281.    return (Result>-EPSILON?(Result<EPSILON?0.0:Result):
  282.    Result);
  283. }
  284.  
  285. /* Normal to a polynomial */
  286. static void normalp(Result, Order, Coeffs, Intersection_Point)
  287. VECTOR *Result;
  288. int Order;
  289. DBL *Coeffs;
  290. VECTOR *Intersection_Point;
  291. {
  292.    int i, xp, yp, zp, wp;
  293.    DBL *a,x[MAX_ORDER+1],y[MAX_ORDER+1],z[MAX_ORDER+1];
  294.    x[0] = 1.0; y[0] = 1.0; z[0] = 1.0;
  295.    x[1] = Intersection_Point->x;
  296.    y[1] = Intersection_Point->y;
  297.    z[1] = Intersection_Point->z;
  298.    for (i=2;i<=Order;i++) {
  299.       x[i] = Intersection_Point->x * x[i-1];
  300.       y[i] = Intersection_Point->y * y[i-1];
  301.       z[i] = Intersection_Point->z * z[i-1];
  302.    }
  303.    a = Coeffs;
  304.    Result->x = 0.0; Result->y = 0.0; Result->z = 0.0;
  305.    for (i=0;i<term_counts[Order];i++) {
  306.       unroll(Order, i, &xp, &yp, &zp, &wp);
  307.       if (xp >= 1) Result->x += xp * a[i] * x[xp-1] * y[yp] * z[zp];
  308.       if (yp >= 1) Result->y += yp * a[i] * x[xp] * y[yp-1] * z[zp];
  309.       if (zp >= 1) Result->z += zp * a[i] * x[xp] * y[yp] * z[zp-1];
  310.    }
  311.    VTemp = sqrt(Result->x * Result->x +
  312.       Result->y * Result->y +
  313.       Result->z * Result->z);
  314.    if (VTemp > 0.0) {
  315.       Result->x /= VTemp;
  316.       Result->y /= VTemp;
  317.       Result->z /= VTemp;
  318.    }
  319.    else {
  320.       Result->x = 1.0;
  321.       Result->y = 0.0;
  322.       Result->z = 0.0;
  323.    }
  324. }
  325.  
  326. static DBL
  327. do_partial_term(q, row, pwr, i, j, k, l)
  328. MATRIX *q;
  329. int row, pwr, i, j, k, l;
  330. {
  331.    DBL result;
  332.    int n;
  333.  
  334.    result = (DBL)(factorials[pwr] /
  335.       (factorials[i]*factorials[j]*factorials[k]*factorials[l]));
  336.    if (i>0) for (n=0;n<i;n++) result *= (*q)[0][row];
  337.    if (j>0) for (n=0;n<j;n++) result *= (*q)[1][row];
  338.    if (k>0) for (n=0;n<k;n++) result *= (*q)[2][row];
  339.    if (l>0) for (n=0;n<l;n++) result *= (*q)[3][row];
  340.    return result;
  341. }
  342.  
  343. /* Using the transformation matrix q, transform the general polynomial
  344.    equation given by a. */
  345.  
  346. static void transform(order, Coeffs, q)
  347. int order;
  348. DBL *Coeffs;
  349. MATRIX *q;
  350. {
  351.    int term_index, partial_index;
  352.    int ip, i, i0, i1, i2, i3;
  353.    int jp, j, j0, j1, j2, j3;
  354.    int kp, k, k0, k1, k2, k3;
  355.    int wp;
  356.    DBL *b, partial_term;
  357.    DBL tempx, tempy, tempz;
  358.  
  359.    for (i=0;i<4;i++) for (j=0;j<4;j++)
  360.       if ((*q)[i][j] > -COEFF_LIMIT && (*q)[i][j] < COEFF_LIMIT)
  361.          (*q)[i][j] = 0.0;
  362.    if ((b = (DBL *)malloc(term_counts[order] * sizeof(DBL))) == NULL) {
  363.       printf("Cannot allocate memory for b in poly transform()\n");
  364.       exit(1);
  365.    }
  366.    for (i=0;i<term_counts[order];i++)
  367.       b[i] = 0.0;
  368.    for (term_index=0;term_index<term_counts[order];term_index++) {
  369.       if (Coeffs[term_index] != 0.0) {
  370.          unroll(order, term_index, &ip, &jp, &kp, &wp);
  371.          /* Step through terms in: (q[0][0]*x+q[0][1]*y+q[0][2]*z+q[0][3])^i */
  372.          for (i=0;i<term_counts[ip];i++) {
  373.             unroll(ip, i, &i0, &i1, &i2, &i3);
  374.             tempx = do_partial_term(q, 0, ip, i0, i1, i2, i3);
  375.             if (tempx != 0.0) {
  376.  
  377.                /* Step through terms in:
  378.             (q[1][0]*x+q[1][1]*y+q[1][2]*z+q[1][3])^j */
  379.                for (j=0;j<term_counts[jp];j++) {
  380.                   unroll(jp, j, &j0, &j1, &j2, &j3);
  381.                   tempy = do_partial_term(q, 1, jp, j0, j1, j2, j3);
  382.                   if (tempy != 0.0) {
  383.  
  384.                      /* Step through terms in:
  385.                   (q[2][0]*x+q[2][1]*y+q[2][2]*z+q[2][3])^k */
  386.                      for (k=0;k<term_counts[kp];k++) {
  387.                         unroll(kp, k, &k0, &k1, &k2, &k3);
  388.                         tempz = do_partial_term(q, 2, kp, k0, k1, k2, k3);
  389.                         if (tempz != 0.0) {
  390.                            /* Figure out it's index, and add into result */
  391.                            partial_index=roll(order,i0+j0+k0,i1+j1+k1,i2+j2+k2);
  392.                            partial_term=Coeffs[term_index]*tempx*tempy*tempz;
  393.                            b[partial_index] += partial_term;
  394.                         }
  395.                      }
  396.                   }
  397.                }
  398.             }
  399.          }
  400.       }
  401.    }
  402.    for (i=0;i<term_counts[order];i++)
  403.       if (b[i] > -1.0e-4 && b[i] < 1.0e-4)
  404.          Coeffs[i] = 0.0;
  405.       else
  406.          Coeffs[i] = b[i];
  407.    free(b);
  408.    /*show_poly("New poly: ", order, Coeffs);
  409. printf("\n");
  410. */
  411. }
  412.  
  413. /* Intersection of a ray and a quartic */
  414. static int intersect_quartic(Ray, Shape, Depths)
  415. RAY *Ray;
  416. POLY *Shape;
  417. DBL *Depths;
  418. {
  419.    DBL x,y,z,x2,y2,z2,x3,y3,z3,x4,y4,z4;
  420.    DBL xx,yy,zz,xx2,yy2,zz2,xx3,yy3,zz3,xx4,yy4,zz4;
  421.    DBL *a, t[5];
  422.    DBL x_z,x_zz,xx_z,xx_zz,x_y,x_yy,xx_y,xx_yy,y_z,y_zz,yy_z,yy_zz,temp;
  423.  
  424.    x = Ray->Initial.x;
  425.    y = Ray->Initial.y;
  426.    z = Ray->Initial.z;
  427.    xx = Ray->Direction.x;
  428.    yy = Ray->Direction.y;
  429.    zz = Ray->Direction.z;
  430.    x2 = x * x;  y2 = y * y;  z2 = z * z;
  431.    x3 = x * x2; y3 = y * y2; z3 = z * z2;
  432.    x4 = x * x3; y4 = y * y3; z4 = z * z3;
  433.    xx2 = xx * xx;  yy2 = yy * yy;  zz2 = zz * zz;
  434.    xx3 = xx * xx2; yy3 = yy * yy2; zz3 = zz * zz2;
  435.    xx4 = xx * xx3; yy4 = yy * yy3; zz4 = zz * zz3;
  436.    a = Shape->Coeffs;
  437.    x_z = x * z; x_zz = x * zz; xx_z = xx * z; xx_zz = xx * zz;
  438.    x_y = x * y; x_yy = x * yy; xx_y = xx * y; xx_yy = xx * yy;
  439.    y_z = y * z; y_zz = y * zz; yy_z = yy * z; yy_zz = yy * zz;
  440.  
  441.    /*
  442.       Determine the coefficients of t^n, where the line is represented
  443.       as (x,y,z) + (xx,yy,zz)*t.
  444.    */
  445.    temp =  a[ 0]*xx4;
  446.    temp += a[ 1]*xx3*yy;
  447.    temp += a[ 2]*xx3*zz;
  448.    temp += a[ 4]*xx2*yy2;
  449.    temp += a[ 5]*xx2*yy_zz;
  450.    temp += a[ 7]*xx2*zz2;
  451.    temp += a[10]*xx*yy3;
  452.    temp += a[11]*xx_zz*yy2;
  453.    temp += a[13]*xx_yy*zz2;
  454.    temp += a[16]*xx*zz3;
  455.    temp += a[20]*yy4;
  456.    temp += a[21]*yy3*zz;
  457.    temp += a[23]*yy2*zz2;
  458.    temp += a[26]*yy*zz3;
  459.    temp += a[30]*zz4;
  460.  
  461.    t[0] = temp;
  462.  
  463.    temp = 4*a[ 0]*x*xx3;
  464.    temp += a[ 1]*(3*xx2*x_yy+xx3*y);
  465.    temp += a[ 2]*(3*xx2*x_zz+xx3*z);
  466.    temp += a[ 3]*xx3;
  467.    temp += a[ 4]*(2*x*xx*yy2+2*xx2*y*yy);
  468.    temp += a[ 5]*(xx2*(y_zz+yy_z)+2*x_yy*xx_zz);
  469.    temp += a[ 6]*xx2*yy;
  470.    temp += a[ 7]*(2*x*xx*zz2+2*xx2*z*zz);
  471.    temp += a[ 8]*xx2*zz;
  472.    temp += a[10]*(x*yy3+3*xx_y*yy2);
  473.    temp += a[11]*(xx*(2*y*yy_zz+yy2*z)+x_zz*yy2);
  474.    temp += a[12]*xx*yy2;
  475.    temp += a[13]*(xx*(y*zz2+2*yy_z*zz)+x_yy*zz2);
  476.    temp += a[14]*xx_yy*zz;
  477.    temp += a[16]*(x*zz3+3*xx_z*zz2);
  478.    temp += a[17]*xx*zz2;
  479.    temp += 4*a[20]*y*yy3;
  480.    temp += a[21]*(3*yy2*y_zz+yy3*z);
  481.    temp += a[22]*yy3;
  482.    temp += zz*(2*a[23]*yy*(y_zz+yy_z)+
  483.       a[24]*yy2+
  484.       zz*(a[26]*(y_zz+3*yy_z)+
  485.       a[27]*yy+zz*(4*a[30]*z+a[31])));
  486.    t[1] = temp;
  487.  
  488.    temp = 6*a[ 0]*x2*xx2;
  489.    temp += 3*a[ 1]*x*xx*(x_yy+xx_y);
  490.    temp += 3*a[ 2]*x*xx*(x_zz+xx_z);
  491.    temp += 3*a[ 3]*x*xx2;
  492.    temp += a[ 4]*(x2*yy2+4*x_yy*xx_y+xx2*y2);
  493.    temp += a[ 5]*(x2*yy_zz+2*x*xx*(y_zz+yy_z)+xx2*y_z);
  494.    temp += a[ 6]*(2*x*xx_yy+xx2*y);
  495.    temp += a[ 7]*(x2*zz2+4*x_zz*xx_z+xx2*z2);
  496.    temp += a[ 8]*(2*x*xx_zz+xx2*z);
  497.    temp += a[ 9]*xx2;
  498.    temp += a[10]*(3*x_y*yy2+3*xx_yy*y2);
  499.    temp += a[11]*(x_yy*(2*y_zz+yy_z)+xx*(y2*zz+2*y*yy_z));
  500.    temp += a[12]*(x*yy2+2*xx_y*yy);
  501.    temp += a[13]*(x_zz*(y_zz+2*yy_z)+xx*(2*y_z*zz+yy*z2));
  502.    temp += a[14]*(x_yy*zz+xx*(y_zz+yy_z));
  503.    temp += a[15]*xx_yy;
  504.    temp += a[16]*(3*x_z*zz2+3*xx_zz*z2);
  505.    temp += a[17]*(x*zz2+2*xx_z*zz);
  506.    temp += a[18]*xx_zz;
  507.    temp += 6*a[20]*y2*yy2;
  508.    temp += 3*a[21]*y*yy*(y_zz+yy_z);
  509.    temp += 3*a[22]*y*yy2;
  510.    temp += a[23]*(y2*zz2+4*y_zz*yy_z+yy2*z2);
  511.    temp += a[24]*(2*y*yy_zz+yy2*z);
  512.    temp += a[25]*yy2;
  513.    temp += zz*(3*a[26]*z*(y_zz+yy_z)+
  514.       a[27]*(y_zz+2*yy_z)+
  515.       a[28]*yy+
  516.       6*a[30]*z2*zz+
  517.       3*a[31]*z*zz+
  518.       a[32]*zz);
  519.    t[2] = temp;
  520.  
  521.    temp = 4*a[ 0]*x3*xx;
  522.    temp += a[ 1]*x2*(x_yy+3*xx_y);
  523.    temp += a[ 2]*x2*(x_zz+3*xx_z);
  524.    temp += 3*a[ 3]*x2*xx;
  525.    temp += 2*a[ 4]*x_y*(x_yy+xx_y);
  526.    temp += a[ 5]*x*(x*(y_zz+yy_z)+2*xx_y*z);
  527.    temp += a[ 6]*x*(x_yy+2*xx_y);
  528.    temp += 2*a[ 7]*x_z*(x_zz+xx_z);
  529.    temp += a[ 8]*x*(x_zz+2*xx_z);
  530.    temp += 2*a[ 9]*x*xx;
  531.    temp += a[10]*(3*x_yy*y2-xx*y3);
  532.    temp += a[11]*(x_y*(y_zz+2*yy_z)+xx_z*y2);
  533.    temp += a[12]*(2*x_y*yy+xx*y2);
  534.    temp += a[13]*(x_z*(2*y_zz+yy_z)+xx_y*z2);
  535.    temp += a[14]*(x*(y_zz+yy_z)+xx_y*z);
  536.    temp += a[15]*(x_yy+xx_y);
  537.    temp += a[16]*(3*x_zz*z2+xx*z3);
  538.    temp += a[17]*(2*x_z*zz+xx*z2);
  539.    temp += a[18]*(x_zz+xx_z);
  540.    temp += a[19]*xx;
  541.    temp += 4*a[20]*y3*yy;
  542.    temp += a[21]*y2*(y_zz+3*yy_z);
  543.    temp += 3*a[22]*y2*yy;
  544.    temp += 2*a[23]*y_z*(y_zz+yy_z);
  545.    temp += a[24]*y*(y_zz+2*yy_z);
  546.    temp += 2*a[25]*y*yy;
  547.    temp += a[26]*(3*y_zz*z2+yy*z3);
  548.    temp += a[27]*(2*y_z*zz+yy*z2);
  549.    temp += a[28]*(y_zz+yy_z);
  550.    temp += a[29]*yy;
  551.    temp += zz*(4*a[30]*z3+3*a[31]*z2+2*a[32]*z+a[33]);
  552.    t[3] = temp;
  553.  
  554.    temp = a[ 0]*x4;
  555.    temp += a[ 1]*x3*y;
  556.    temp += a[ 2]*x3*z;
  557.    temp += a[ 3]*x3;
  558.    temp += a[ 4]*x2*y2;
  559.    temp += a[ 5]*x2*y_z;
  560.    temp += a[ 6]*x2*y;
  561.    temp += a[ 7]*x2*z2;
  562.    temp += a[ 8]*x2*z;
  563.    temp += a[ 9]*x2;
  564.    temp += a[10]*x*y3;
  565.    temp += a[11]*x*y2*z;
  566.    temp += a[12]*x*y2;
  567.    temp += a[13]*x_y*z2;
  568.    temp += a[14]*x_y*z;
  569.    temp += a[15]*x_y;
  570.    temp += a[16]*x*z3;
  571.    temp += a[17]*x*z2;
  572.    temp += a[18]*x_z;
  573.    temp += a[19]*x;
  574.    temp += a[20]*y4;
  575.    temp += a[21]*y3*z;
  576.    temp += a[22]*y3;
  577.    temp += a[23]*y2*z2;
  578.    temp += a[24]*y2*z;
  579.    temp += a[25]*y2;
  580.    temp += a[26]*y*z3;
  581.    temp += a[27]*y*z2;
  582.    temp += a[28]*y_z;
  583.    temp += a[29]*y;
  584.    temp += a[30]*z4;
  585.    temp += a[31]*z3;
  586.    temp += a[32]*z2;
  587.    temp += a[33]*z;
  588.    temp += a[34];
  589.    t[4] = temp;
  590.  
  591.    if (Shape->Sturm_Flag)
  592.       if (t[0] == 0.0)
  593.          if (t[1] == 0.0)
  594.             return solve_quadratic(&t[2], Depths);
  595.          else
  596.             return polysolve(3, &t[1], Depths);
  597.       else
  598.          return polysolve(4, &t[0], Depths);
  599.    else
  600.       return solve_quartic(&t[0], Depths);
  601. }
  602.  
  603. #ifdef TEST
  604. /* Normal to a quartic */
  605. static void quartic_normal(Result, Object, Intersection_Point)
  606. VECTOR *Result, *Intersection_Point;
  607. OBJECT *Object;
  608. {
  609.    POLY *Shape = (POLY *) Object;
  610.    DBL *a,x,y,z,x2,y2,z2,x3,y3,z3,temp;
  611.  
  612.    a = Shape->Coeffs;
  613.    x = Intersection_Point->x;
  614.    y = Intersection_Point->y;
  615.    z = Intersection_Point->z;
  616.    x2 = x * x;  y2 = y * y;  z2 = z * z;
  617.    x3 = x * x2; y3 = y * y2; z3 = z * z2;
  618.  
  619.    temp = a[ 4]*y2+y*(a[ 5]*z+a[ 6])+a[ 7]*z2+a[ 8]*z+a[ 9];
  620.    Result->x = 4*a[ 0]*x3+3*x2*(a[ 1]*y+a[ 2]*z+a[ 3])+
  621.    2*x*(temp)+
  622.    a[10]*y3+y2*(a[11]*z+a[12])+y*(a[13]*z2+a[14]*z+a[15])+
  623.    a[16]*z3+a[17]*z2+a[18]*z+a[19];
  624.  
  625.    temp = 3*a[10]*y2+2*y*(a[11]*z+a[12])+a[13]*z2+a[14]*z+a[15];
  626.    Result->y = a[ 1]*x3+x2*(2*a[ 4]*y+a[ 5]*z+a[ 6])+
  627.    x*(temp)+
  628.    4*a[20]*y3+3*y2*(a[21]*z+a[22])+2*y*(a[23]*z2+a[24]*z+a[25])+
  629.    a[26]*z3+a[27]*z2+a[28]*z+a[29];
  630.  
  631.    temp = a[11]*y2+y*(2*a[13]*z+a[14])+3*a[16]*z2+2*a[17]*z+a[18];
  632.    Result->z = a[ 2]*x3+x2*(a[ 5]*y+2*a[ 7]*z+a[ 8])+
  633.    x*(temp)+
  634.    a[21]*y3+y2*(2*a[23]*z+a[24])+y*(3*a[26]*z2+2*a[27]*z+a[28])+
  635.    4*a[30]*z3+3*a[31]*z2+2*a[32]*z+a[33];
  636.  
  637.    VTemp=Result->x * Result->x + Result->y * Result->y + Result->z * Result->z;
  638.    VTemp=sqrt(VTemp);
  639.    if (VTemp > 0.0) {
  640.       Result->x /= VTemp;
  641.       Result->y /= VTemp;
  642.       Result->z /= VTemp;
  643.    }
  644.    else {
  645.       Result->x = 1.0;
  646.       Result->y = 0.0;
  647.       Result->z = 0.0;
  648.    }
  649. }
  650. #endif
  651.  
  652. /* Normal to a quartic */
  653. static void quartic_normal(Result, Object, Intersection_Point)
  654. VECTOR *Result, *Intersection_Point;
  655. OBJECT *Object;
  656. {
  657.    POLY *Shape = (POLY *) Object;
  658.    DBL *a,x,y,z,x2,y2,z2,x3,y3,z3;
  659.  
  660.    a = Shape->Coeffs;
  661.    x = Intersection_Point->x;
  662.    y = Intersection_Point->y;
  663.    z = Intersection_Point->z;
  664.    x2 = x * x;  y2 = y * y;  z2 = z * z;
  665.    x3 = x * x2; y3 = y * y2; z3 = z * z2;
  666.  
  667.    Result->x = 4*a[ 0]*x3+3*x2*(a[ 1]*y+a[ 2]*z+a[ 3])+
  668.    2*x*(a[ 4]*y2+y*(a[ 5]*z+a[ 6])+a[ 7]*z2+a[ 8]*z+a[ 9])+
  669.    a[10]*y3+y2*(a[11]*z+a[12])+y*(a[13]*z2+a[14]*z+a[15])+
  670.    a[16]*z3+a[17]*z2+a[18]*z+a[19];
  671.  
  672.    Result->y = a[ 1]*x3+x2*(2*a[ 4]*y+a[ 5]*z+a[ 6])+
  673.    x*(3*a[10]*y2+2*y*(a[11]*z+a[12])+a[13]*z2+a[14]*z+a[15])+
  674.    4*a[20]*y3+3*y2*(a[21]*z+a[22])+2*y*(a[23]*z2+a[24]*z+a[25])+
  675.    a[26]*z3+a[27]*z2+a[28]*z+a[29];
  676.  
  677.    Result->z = a[ 2]*x3+x2*(a[ 5]*y+2*a[ 7]*z+a[ 8])+
  678.    x*(a[11]*y2+y*(2*a[13]*z+a[14])+3*a[16]*z2+2*a[17]*z+a[18])+
  679.    a[21]*y3+y2*(2*a[23]*z+a[24])+y*(3*a[26]*z2+2*a[27]*z+a[28])+
  680.    4*a[30]*z3+3*a[31]*z2+2*a[32]*z+a[33];
  681.    VTemp = sqrt(Result->x*Result->x+Result->y*Result->y+Result->z*Result->z);
  682.    if (VTemp > 0.0) {
  683.       Result->x /= VTemp;
  684.       Result->y /= VTemp;
  685.       Result->z /= VTemp;
  686.    }
  687.    else {
  688.       Result->x = 1.0;
  689.       Result->y = 0.0;
  690.       Result->z = 0.0;
  691.    }
  692. }
  693.  
  694. int Inside_Poly (Test_Point, Object)
  695. VECTOR *Test_Point;
  696. OBJECT *Object;
  697. {
  698.    VECTOR New_Point;
  699.    POLY *Shape = (POLY *) Object;
  700.    DBL Result;
  701.  
  702.    /* Transform the point into polynomial's space */
  703.    if (Shape->Transform != NULL)
  704.       MInverseTransformVector(&New_Point, Test_Point, Shape->Transform);
  705.    else
  706.       New_Point = *Test_Point;
  707.  
  708.    Result = inside(&New_Point, Shape->Order, Shape->Coeffs);
  709.    if (Result < Small_Tolerance)
  710.       return ((int)(1-Shape->Inverted));
  711.    else
  712.       return ((int)Shape->Inverted);
  713. }
  714.  
  715. /* Normal to a polynomial */
  716. void Poly_Normal(Result, Object, Intersection_Point)
  717. VECTOR *Result, *Intersection_Point;
  718. OBJECT *Object;
  719. {
  720.    POLY *Shape = (POLY *) Object;
  721.    VECTOR New_Point;
  722.  
  723.    /* Transform the point into the polynomials space */
  724.    if (Shape->Transform != NULL)
  725.       MInverseTransformVector(&New_Point, Intersection_Point, Shape->Transform);
  726.    else {
  727.       New_Point.x = Intersection_Point->x;
  728.       New_Point.y = Intersection_Point->y;
  729.       New_Point.z = Intersection_Point->z;
  730.    }
  731.  
  732.       if (Shape->Order == 4)
  733.          quartic_normal(Result, Object, &New_Point);
  734.       else
  735.          normalp(Result, Shape->Order, Shape->Coeffs, &New_Point);
  736.  
  737.    /* Transform back to world space */
  738.    if (Shape->Transform != NULL)
  739.       MTransNormal(Result, Result, Shape->Transform);
  740.    VNormalize(*Result, *Result);
  741. }
  742.  
  743. /* Make a copy of a polynomial object */
  744. void *Copy_Poly(Object)
  745. OBJECT *Object;
  746. {
  747.    POLY *Shape = (POLY *)Object;
  748.    POLY *New_Shape = Get_Poly_Shape(Shape->Order);
  749.    int i;
  750.  
  751.    New_Shape->Shape_Texture = Shape->Shape_Texture;
  752.    New_Shape->Shape_Colour = Shape->Shape_Colour;
  753.    New_Shape->Next_Object = NULL;
  754.    New_Shape->Sturm_Flag = Shape->Sturm_Flag;
  755.    New_Shape->Inverted = Shape->Inverted;
  756.  
  757.    /* Copy any associated transformation */
  758.    if (Shape->Transform != NULL) {
  759.       New_Shape->Transform = Get_Transformation();
  760.       memcpy(New_Shape->Transform, Shape->Transform, sizeof(TRANSFORMATION));
  761.    }
  762.    for (i=0;i<term_counts[New_Shape->Order];i++)
  763.       New_Shape->Coeffs[i] = Shape->Coeffs[i];
  764.  
  765.    /* Copy any associated texture */
  766.    if (Shape->Shape_Texture != NULL)
  767.       New_Shape->Shape_Texture = Copy_Texture (Shape->Shape_Texture);
  768.  
  769.    return (void *)(New_Shape);
  770. }
  771.  
  772. void Translate_Poly (Object, Vector)
  773. OBJECT *Object;
  774. VECTOR *Vector;
  775. {
  776.    TRANSFORMATION Transform;
  777.    POLY *Shape = (POLY *)Object;
  778.    if (Shape->Transform == NULL)
  779.       Shape->Transform = Get_Transformation();
  780.    Get_Translation_Transformation(&Transform, Vector);
  781.    Compose_Transformations(Shape->Transform, &Transform);
  782.  
  783.    Translate_Texture (&Shape->Shape_Texture, Vector);
  784. }
  785.  
  786. void Rotate_Poly (Object, Vector)
  787. OBJECT *Object;
  788. VECTOR *Vector;
  789. {
  790.    TRANSFORMATION Transform;
  791.    POLY *Shape = (POLY *)Object;
  792.    if (Shape->Transform == NULL)
  793.       Shape->Transform = Get_Transformation();
  794.    Get_Rotation_Transformation(&Transform, Vector);
  795.    Compose_Transformations(Shape->Transform, &Transform);
  796.  
  797.    Rotate_Texture (&Shape->Shape_Texture, Vector);
  798. }
  799.  
  800. void Scale_Poly (Object, Vector)
  801. OBJECT *Object;
  802. VECTOR *Vector;
  803. {
  804.    TRANSFORMATION Transform;
  805.    POLY *Shape = (POLY *)Object;
  806.    if (Shape->Transform == NULL)
  807.       Shape->Transform = Get_Transformation();
  808.    Get_Scaling_Transformation(&Transform, Vector);
  809.    Compose_Transformations(Shape->Transform, &Transform);
  810.  
  811.    Scale_Texture (&Shape->Shape_Texture, Vector);
  812. }
  813.  
  814. void Invert_Poly (Object)
  815. OBJECT *Object;
  816. {
  817.    ((POLY *) Object)->Inverted = 1 - ((POLY *)Object)->Inverted;
  818. }
  819.