home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / mac / raytrace / dkb / mcprt102.bnh / DKBTrace-Mac / quartics.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-07-04  |  17.4 KB  |  557 lines
  1. /*****************************************************************************
  2. *
  3. *                                   quartics.c
  4. *
  5. *   from DKBTrace (c) 1990  David Buck
  6. *
  7. *  This module implements the code for the quartic (4th order) shapes.
  8. *
  9. *  This file was written by Alexander Enzmann.  He wrote the code for DKB 2.11
  10. *  QUARTICS (4th order shapes) and generously provided us these enhancements.
  11. *
  12. * This software is freely distributable. The source and/or object code may be
  13. * copied or uploaded to communications services so long as this notice remains
  14. * at the top of each file.  If any changes are made to the program, you must
  15. * clearly indicate in the documentation and in the programs startup message
  16. * who it was who made the changes. The documentation should also describe what
  17. * those changes were. This software may not be included in whole or in
  18. * part into any commercial package without the express written consent of the
  19. * author.  It may, however, be included in other public domain or freely
  20. * distributed software so long as the proper credit for the software is given.
  21. *
  22. * This software is provided as is without any guarantees or warranty. Although
  23. * the author has attempted to find and correct any bugs in the software, he
  24. * is not responsible for any damage caused by the use of the software.  The
  25. * author is under no obligation to provide service, corrections, or upgrades
  26. * to this package.
  27. *
  28. * Despite all the legal stuff above, if you do find bugs, I would like to hear
  29. * about them.  Also, if you have any comments or questions, you may contact me
  30. * at the following address:
  31. *
  32. *     David Buck
  33. *     22C Sonnet Cres.
  34. *     Nepean Ontario
  35. *     Canada, K2H 8W7
  36. *
  37. *  I can also be reached on the following bulleton boards:
  38. *
  39. *     OMX              (613) 731-3419
  40. *     Mystic           (613) 596-4249  or  (613) 596-4772
  41. *
  42. *  Fidonet:   1:163/109.9
  43. *  Internet:  dbuck@ccs.carleton.ca
  44. *  The "You Can Call Me RAY" BBS    (708) 358-5611
  45. *
  46. *  IBM Port by Aaron A. Collins. Aaron may be reached on the following BBS'es:
  47. *
  48. *     The "You Can Call Me RAY" BBS (708) 358-5611
  49. *     The Information Exchange BBS  (708) 945-5575
  50. *
  51. *****************************************************************************/
  52.  
  53. #include "frame.h"
  54. #include "vector.h"
  55. #include "dkbproto.h"
  56.  
  57.  
  58. /* Basic form of a quartic equation
  59.    a00*x^4 + a01*x^3*y + a02*x^3*z + a03*x^3 + a04*x^2*y^2 + 
  60.    a05*x^2*y*z + a06*x^2*y + a07*x^2*z^2 + a08*x^2*z + a09*x^2 + 
  61.    a10*x*y^3 + a11*x*y^2*z + a12*x*y^2 + a13*x*y*z^2 + a14*x*y*z + 
  62.    a15*x*y + a16*x*z^3 + a17*x*z^2 + a18*x*z + a19*x + a20*y^4 + 
  63.    a21*y^3*z + a22*y^3 + a23*y^2*z^2 + a24*y^2*z + a25*y^2 + a26*y*z^3 + 
  64.    a27*y*z^2 + a28*y*z + a29*y + a30*z^4 + a31*z^3 + a32*z^2 + a33*z + a34
  65. */
  66.  
  67. int factorials[5] = {1, 1, 2, 6, 24};    /* 0!, 1!, 2!, 3!, 4! */
  68. int term_counts[5] = {1, 4, 10, 20, 35};
  69. int terms4[35][4] = {
  70.    {4, 0, 0, 0 }, {3, 1, 0, 0 }, {3, 0, 1, 0 }, {3, 0, 0, 1 }, {2, 2, 0, 0 },
  71.    {2, 1, 1, 0 }, {2, 1, 0, 1 }, {2, 0, 2, 0 }, {2, 0, 1, 1 }, {2, 0, 0, 2 },
  72.    {1, 3, 0, 0 }, {1, 2, 1, 0 }, {1, 2, 0, 1 }, {1, 1, 2, 0 }, {1, 1, 1, 1 },
  73.    {1, 1, 0, 2 }, {1, 0, 3, 0 }, {1, 0, 2, 1 }, {1, 0, 1, 2 }, {1, 0, 0, 3 },
  74.    {0, 4, 0, 0 }, {0, 3, 1, 0 }, {0, 3, 0, 1 }, {0, 2, 2, 0 }, {0, 2, 1, 1 },
  75.    {0, 2, 0, 2 }, {0, 1, 3, 0 }, {0, 1, 2, 1 }, {0, 1, 1, 2 }, {0, 1, 0, 3 },
  76.    {0, 0, 4, 0 }, {0, 0, 3, 1 }, {0, 0, 2, 2 }, {0, 0, 1, 3 }, {0, 0, 0, 4 },
  77.    };
  78.  
  79. int terms3[20][4] = {
  80.    {3, 0, 0, 0 }, {2, 1, 0, 0 }, {2, 0, 1, 0 }, {2, 0, 0, 1 }, {1, 2, 0, 0 },
  81.    {1, 1, 1, 0 }, {1, 1, 0, 1 }, {1, 0, 2, 0 }, {1, 0, 1, 1 }, {1, 0, 0, 2 },
  82.    {0, 3, 0, 0 }, {0, 2, 1, 0 }, {0, 2, 0, 1 }, {0, 1, 2, 0 }, {0, 1, 1, 1 },
  83.    {0, 1, 0, 2 }, {0, 0, 3, 0 }, {0, 0, 2, 1 }, {0, 0, 1, 2 }, {0, 0, 0, 3 },
  84.    };
  85.  
  86. int terms2[10][4] = {
  87.    {2, 0, 0, 0 }, {1, 1, 0, 0 }, {1, 0, 1, 0 }, {1, 0, 0, 1 }, {0, 2, 0, 0 },
  88.    {0, 1, 1, 0 }, {0, 1, 0, 1 }, {0, 0, 2, 0 }, {0, 0, 1, 1 }, {0, 0, 0, 2 },
  89.    };
  90.  
  91. int terms1[4][4] = {
  92.    {1, 0, 0, 0 }, {0, 1, 0, 0 }, {0, 0, 1, 0 }, {0, 0, 0, 1 },
  93.    };
  94.  
  95.  
  96. METHODS Quartic_Methods =
  97.    { Object_Intersect, All_Quartic_Intersections,
  98.      Inside_Quartic, Quartic_Normal, Copy_Quartic,
  99.      Translate_Quartic, Rotate_Quartic,
  100.      Scale_Quartic, Invert_Quartic};
  101.  
  102. extern long Ray_Quartic_Tests, Ray_Quartic_Tests_Succeeded;
  103.  
  104. int
  105. All_Quartic_Intersections(Object, Ray, Depth_Queue)
  106.    OBJECT *Object;
  107.    RAY *Ray;
  108.    PRIOQ *Depth_Queue;
  109. {
  110.    QUARTIC *Shape = (QUARTIC *) Object;
  111.    DBL Depths[4];
  112.    VECTOR Intersection_Point;
  113.    INTERSECTION Local_Element;
  114.    int cnt, i, j, Intersection_Found;
  115.  
  116.    Intersection_Found = FALSE;
  117.    Ray_Quartic_Tests++;
  118.    cnt = Intersect_Quartic(Ray, Shape, &Depths[0]);
  119.    if (cnt > 0) Ray_Quartic_Tests_Succeeded++;
  120.    for (i=0;i<cnt;i++) {
  121.       for (j=0;j<i;j++)
  122.      if (Depths[i] == Depths[j]) goto l0;
  123.  
  124.       if (Depths[i] < Small_Tolerance)
  125.          continue;
  126.  
  127.       Local_Element.Depth = Depths[i];
  128.       Local_Element.Object = Shape -> Parent_Object;
  129.       VScale (Intersection_Point, Ray -> Direction, Depths[i]);
  130.       VAdd (Intersection_Point, Intersection_Point, Ray -> Initial);
  131.       Local_Element.Point = Intersection_Point;
  132.       Local_Element.Shape = (SHAPE *)Shape;
  133.       pq_add (Depth_Queue, &Local_Element);
  134.       Intersection_Found = TRUE;
  135.       l0:;
  136.       }
  137.    return (Intersection_Found);
  138. }
  139.  
  140. /* Intersection of a ray and a quartic */
  141. int
  142. Intersect_Quartic(Ray, Shape, Depths)
  143.    RAY *Ray;
  144.    QUARTIC *Shape;
  145.    DBL *Depths;
  146. {
  147.    DBL x,y,z,x2,y2,z2,x3,y3,z3,x4,y4,z4;
  148.    DBL xx,yy,zz,xx2,yy2,zz2,xx3,yy3,zz3,xx4,yy4,zz4;
  149.    DBL *a, t[5];
  150.    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;
  151.  
  152.    x = Ray->Initial.x;
  153.    y = Ray->Initial.y;
  154.    z = Ray->Initial.z;
  155.    xx = Ray->Direction.x;
  156.    yy = Ray->Direction.y;
  157.    zz = Ray->Direction.z;
  158.    x2 = x * x;  y2 = y * y;  z2 = z * z;
  159.    x3 = x * x2; y3 = y * y2; z3 = z * z2;
  160.    x4 = x * x3; y4 = y * y3; z4 = z * z3;
  161.    xx2 = xx * xx;  yy2 = yy * yy;  zz2 = zz * zz;
  162.    xx3 = xx * xx2; yy3 = yy * yy2; zz3 = zz * zz2;
  163.    xx4 = xx * xx3; yy4 = yy * yy3; zz4 = zz * zz3;
  164.    a = Shape->Coeffs;
  165.    x_z = x * z; x_zz = x * zz; xx_z = xx * z; xx_zz = xx * zz;
  166.    x_y = x * y; x_yy = x * yy; xx_y = xx * y; xx_yy = xx * yy;
  167.    y_z = y * z; y_zz = y * zz; yy_z = yy * z; yy_zz = yy * zz;
  168.  
  169.    /*
  170.       Determine the coefficients of t^n, where the line is represented
  171.       as (x,y,z) + (xx,yy,zz)*t.
  172.    */
  173.     t[0] = (a[ 0]*xx4+a[ 1]*xx3*yy+a[ 2]*xx3*zz+a[ 4]*xx2*yy2+
  174.             a[ 5]*xx2*yy_zz+a[ 7]*xx2*zz2+a[10]*xx*yy3+
  175.             a[11]*xx_zz*yy2+a[13]*xx_yy*zz2+a[16]*xx*zz3+
  176.             a[20]*yy4+a[21]*yy3*zz+a[23]*yy2*zz2+a[26]*yy*zz3+a[30]*zz4);
  177.  
  178.     temp = 4*a[ 0]*x*xx3;
  179.     temp += a[ 1]*(3*xx2*x_yy+xx3*y);
  180.     temp += a[ 2]*(3*xx2*x_zz+xx3*z);
  181.     temp += a[ 3]*xx3;
  182.     temp += a[ 4]*(2*x*xx*yy2+2*xx2*y*yy);
  183.     temp += a[ 5]*(xx2*(y_zz+yy_z)+2*x_yy*xx_zz);
  184.     temp += a[ 6]*xx2*yy;
  185.     temp += a[ 7]*(2*x*xx*zz2+2*xx2*z*zz);
  186.     temp += a[ 8]*xx2*zz;
  187.     temp += a[10]*(x*yy3+3*xx_y*yy2);
  188.     temp += a[11]*(xx*(2*y*yy_zz+yy2*z)+x_zz*yy2);
  189.     temp += a[12]*xx*yy2;
  190.     temp += a[13]*(xx*(y*zz2+2*yy_z*zz)+x_yy*zz2);
  191.     temp += a[14]*xx_yy*zz;
  192.     temp += a[16]*(x*zz3+3*xx_z*zz2);
  193.     temp += a[17]*xx*zz2;
  194.     temp += 4*a[20]*y*yy3;
  195.     temp += a[21]*(3*yy2*y_zz+yy3*z);
  196.     temp += a[22]*yy3;
  197.     temp += zz*(2*a[23]*yy*(y_zz+yy_z)+
  198.             a[24]*yy2+
  199.             zz*(a[26]*(y_zz+3*yy_z)+
  200.             a[27]*yy+zz*(4*a[30]*z+a[31])));
  201.     t[1] = temp;
  202.  
  203.     temp = 6*a[ 0]*x2*xx2;
  204.     temp += 3*a[ 1]*x*xx*(x_yy+xx_y);
  205.     temp += 3*a[ 2]*x*xx*(x_zz+xx_z);
  206.     temp += 3*a[ 3]*x*xx2;
  207.     temp += a[ 4]*(x2*yy2+4*x_yy*xx_y+xx2*y2);
  208.     temp += a[ 5]*(x2*yy_zz+2*x*xx*(y_zz+yy_z)+xx2*y_z);
  209.     temp += a[ 6]*(2*x*xx_yy+xx2*y);
  210.     temp += a[ 7]*(x2*zz2+4*x_zz*xx_z+xx2*z2);
  211.     temp += a[ 8]*(2*x*xx_zz+xx2*z);
  212.     temp += a[ 9]*xx2;
  213.     temp += a[10]*(3*x_y*yy2+3*xx_yy*y2);
  214.     temp += a[11]*(x_yy*(2*y_zz+yy_z)+xx*(y2*zz+2*y*yy_z));
  215.     temp += a[12]*(x*yy2+2*xx_y*yy);
  216.     temp += a[13]*(x_zz*(y_zz+2*yy_z)+xx*(2*y_z*zz+yy*z2));
  217.     temp += a[14]*(x_yy*zz+xx*(y_zz+yy_z));
  218.     temp += a[15]*xx_yy;
  219.     temp += a[16]*(3*x_z*zz2+3*xx_zz*z2);
  220.     temp += a[17]*(x*zz2+2*xx_z*zz);
  221.     temp += a[18]*xx_zz;
  222.     temp += 6*a[20]*y2*yy2;
  223.     temp += 3*a[21]*y*yy*(y_zz+yy_z);
  224.     temp += 3*a[22]*y*yy2;
  225.     temp += a[23]*(y2*zz2+4*y_zz*yy_z+yy2*z2);
  226.     temp += a[24]*(2*y*yy_zz+yy2*z);
  227.     temp += a[25]*yy2;
  228.     temp += zz*(3*a[26]*z*(y_zz+yy_z)+
  229.             a[27]*(y_zz+2*yy_z)+
  230.             a[28]*yy+
  231.             6*a[30]*z2*zz+
  232.             3*a[31]*z*zz+
  233.             a[32]*zz);
  234.     t[2] = temp;
  235.  
  236.     temp = 4*a[ 0]*x3*xx;
  237.     temp += a[ 1]*x2*(x_yy+3*xx_y);
  238.     temp += a[ 2]*x2*(x_zz+3*xx_z);
  239.     temp += 3*a[ 3]*x2*xx;
  240.     temp += 2*a[ 4]*x_y*(x_yy+xx_y);
  241.     temp += a[ 5]*x*(x*(y_zz+yy_z)+2*xx_y*z);
  242.     temp += a[ 6]*x*(x_yy+2*xx_y);
  243.     temp += 2*a[ 7]*x_z*(x_zz+xx_z);
  244.     temp += a[ 8]*x*(x_zz+2*xx_z);
  245.     temp += 2*a[ 9]*x*xx;
  246.     temp += a[10]*(3*x_yy*y2-xx*y3);
  247.     temp += a[11]*(x_y*(y_zz+2*yy_z)+xx_z*y2);
  248.     temp += a[12]*(2*x_y*yy+xx*y2);
  249.     temp += a[13]*(x_z*(2*y_zz+yy_z)+xx_y*z2);
  250.     temp += a[14]*(x*(y_zz+yy_z)+xx_y*z);
  251.     temp += a[15]*(x_yy+xx_y);
  252.     temp += a[16]*(3*x_zz*z2+xx*z3);
  253.     temp += a[17]*(2*x_z*zz+xx*z2);
  254.     temp += a[18]*(x_zz+xx_z);
  255.     temp += a[19]*xx;
  256.     temp += 4*a[20]*y3*yy;
  257.     temp += a[21]*y2*(y_zz+3*yy_z);
  258.     temp += 3*a[22]*y2*yy;
  259.     temp += 2*a[23]*y_z*(y_zz+yy_z);
  260.     temp += a[24]*y*(y_zz+2*yy_z);
  261.     temp += 2*a[25]*y*yy;
  262.     temp += a[26]*(3*y_zz*z2+yy*z3);
  263.     temp += a[27]*(2*y_z*zz+yy*z2);
  264.     temp += a[28]*(y_zz+yy_z);
  265.     temp += a[29]*yy;
  266.     temp += zz*(4*a[30]*z3+3*a[31]*z2+2*a[32]*z+a[33]);
  267.     t[3] = temp;
  268.  
  269.     t[4] = a[ 0]*x4+a[ 1]*x3*y+a[ 2]*x3*z+a[ 3]*x3+a[ 4]*x2*y2+
  270.             a[ 5]*x2*y_z+a[ 6]*x2*y+a[ 7]*x2*z2+a[ 8]*x2*z+a[ 9]*x2+
  271.             a[10]*x*y3+a[11]*x*y2*z+a[12]*x*y2+a[13]*x_y*z2+a[14]*x_y*z+
  272.             a[15]*x_y+a[16]*x*z3+a[17]*x*z2+a[18]*x_z+a[19]*x+a[20]*y4+
  273.             a[21]*y3*z+a[22]*y3+a[23]*y2*z2+a[24]*y2*z+a[25]*y2+
  274.             a[26]*y*z3+a[27]*y*z2+a[28]*y_z+a[29]*y+a[30]*z4+a[31]*z3+
  275.             a[32]*z2+a[33]*z+a[34];
  276.             
  277.    return solve_quartic(t, Depths);
  278. }
  279.  
  280. int
  281. Inside_Quartic (point, Object)
  282.    VECTOR *point;
  283.    OBJECT *Object;
  284. {
  285.    QUARTIC *Shape = (QUARTIC *) Object;
  286.    DBL Result,*a,x,y,z,x2,y2,z2,x3,y3,z3,x4,y4,z4;
  287.  
  288.    x = point->x; y = point->y; z = point->z;
  289.    x2 = x * x;  y2 = y * y;  z2 = z * z;
  290.    x3 = x * x2; y3 = y * y2; z3 = z * z2;
  291.    x4 = x * x3; y4 = y * y3; z4 = z * z3;
  292.    a = Shape->Coeffs;
  293.  
  294.    Result = a[ 0]*x4+a[ 1]*x3*y+a[ 2]*x3*z+a[ 3]*x3+a[ 4]*x2*y2+
  295.         a[ 5]*x2*y*z+a[ 6]*x2*y+a[ 7]*x2*z2+a[ 8]*x2*z+a[ 9]*x2+
  296.         a[10]*x*y3+a[11]*x*y2*z+a[12]*x*y2+a[13]*x*y*z2+a[14]*x*y*z+
  297.         a[15]*x*y+a[16]*x*z3+a[17]*x*z2+a[18]*x*z+a[19]*x+a[20]*y4+
  298.         a[21]*y3*z+a[22]*y3+a[23]*y2*z2+a[24]*y2*z+a[25]*y2+
  299.         a[26]*y*z3+a[27]*y*z2+a[28]*y*z+a[29]*y+a[30]*z4+a[31]*z3+
  300.         a[32]*z2+a[33]*z+a[34];
  301.  
  302.    if (Result < Small_Tolerance)
  303.       return (TRUE);
  304.    else
  305.       return (FALSE);
  306. }
  307.  
  308. /* Normal to a quartic */
  309. void
  310. Quartic_Normal(Result, Object, Intersection_Point)
  311.    VECTOR *Result, *Intersection_Point;
  312.    OBJECT *Object;
  313. {
  314.    QUARTIC *Shape = (QUARTIC *) Object;
  315.    DBL *a,x,y,z,x2,y2,z2,x3,y3,z3,Len,temp;
  316.  
  317.    a = Shape->Coeffs;
  318.    x = Intersection_Point->x;
  319.    y = Intersection_Point->y;
  320.    z = Intersection_Point->z;
  321.    x2 = x * x;  y2 = y * y;  z2 = z * z;
  322.    x3 = x * x2; y3 = y * y2; z3 = z * z2;
  323.  
  324.    temp = a[ 4]*y2+y*(a[ 5]*z+a[ 6])+a[ 7]*z2+a[ 8]*z+a[ 9];
  325.    Result->x = 4*a[ 0]*x3+3*x2*(a[ 1]*y+a[ 2]*z+a[ 3])+
  326.            2*x*(temp)+
  327.            a[10]*y3+y2*(a[11]*z+a[12])+y*(a[13]*z2+a[14]*z+a[15])+
  328.            a[16]*z3+a[17]*z2+a[18]*z+a[19];
  329.  
  330.    temp = 3*a[10]*y2+2*y*(a[11]*z+a[12])+a[13]*z2+a[14]*z+a[15];
  331.    Result->y = a[ 1]*x3+x2*(2*a[ 4]*y+a[ 5]*z+a[ 6])+
  332.            x*(temp)+
  333.            4*a[20]*y3+3*y2*(a[21]*z+a[22])+2*y*(a[23]*z2+a[24]*z+a[25])+
  334.            a[26]*z3+a[27]*z2+a[28]*z+a[29];
  335.  
  336.    temp = a[11]*y2+y*(2*a[13]*z+a[14])+3*a[16]*z2+2*a[17]*z+a[18];
  337.    Result->z = a[ 2]*x3+x2*(a[ 5]*y+2*a[ 7]*z+a[ 8])+
  338.            x*(temp)+
  339.            a[21]*y3+y2*(2*a[23]*z+a[24])+y*(3*a[26]*z2+2*a[27]*z+a[28])+
  340.            4*a[30]*z3+3*a[31]*z2+2*a[32]*z+a[33];
  341.  
  342.    Len = Result->x * Result->x + Result->y * Result->y + Result->z * Result->z;
  343.    Len = sqrt(Len);
  344.    if (Len == 0.0) {
  345.       /* The normal is not defined at this point of the surface.  Set it
  346.          to any arbitrary direction. */
  347.       Result->x = 1.0;
  348.       Result->y = 0.0;
  349.       Result->z = 0.0;
  350.       }
  351.    else {
  352.       Result->x /= Len;        /* normalize the normal */
  353.       Result->y /= Len;
  354.       Result->z /= Len;
  355.       }
  356.    }
  357.  
  358. /* Make a copy of a quartic object */
  359. void *
  360. Copy_Quartic(Object)
  361.    OBJECT *Object;
  362. {
  363.    QUARTIC *New_Shape;
  364.  
  365.    New_Shape = Get_Quartic_Shape ();
  366.    *New_Shape = *((QUARTIC *) Object);
  367.    New_Shape -> Next_Object = NULL;
  368.  
  369.    if (New_Shape->Shape_Texture != NULL)
  370.       New_Shape->Shape_Texture = Copy_Texture (New_Shape->Shape_Texture);
  371.  
  372.    return (New_Shape);
  373. }
  374.  
  375. /* Move a quartic */
  376. void
  377. Translate_Quartic(Object, Vector)
  378.    OBJECT *Object;
  379.    VECTOR *Vector;
  380. {
  381.    TRANSFORMATION Transformation;
  382.    Get_Translation_Transformation(&Transformation, Vector);
  383.    Transform_Quartic((QUARTIC *)Object,
  384.              (MATRIX *)&(Transformation.inverse[0][0]));
  385.    Translate_Texture (&((QUARTIC *) Object)->Shape_Texture, Vector);
  386. }
  387.  
  388. /* Rotation of a quartic */
  389. void
  390. Rotate_Quartic(Object, Vector)
  391.    OBJECT *Object;
  392.    VECTOR *Vector;
  393. {
  394.    TRANSFORMATION Transformation;
  395.    Get_Rotation_Transformation(&Transformation, Vector);
  396.    Transform_Quartic((QUARTIC *)Object,
  397.              (MATRIX *)&(Transformation.inverse[0][0]));
  398.    Rotate_Texture (&((QUARTIC *) Object)->Shape_Texture, Vector);
  399. }
  400.  
  401. /* Resize a quartic */
  402. void
  403. Scale_Quartic(Object, Vector)
  404.    OBJECT *Object;
  405.    VECTOR *Vector;
  406. {
  407.    TRANSFORMATION Transformation;
  408.    Get_Scaling_Transformation (&Transformation, Vector);
  409.    Transform_Quartic((QUARTIC *)Object,
  410.              (MATRIX *)&(Transformation.inverse[0][0]));
  411.    Scale_Texture (&((QUARTIC *) Object)->Shape_Texture, Vector);
  412. }
  413.  
  414. void
  415. Invert_Quartic(Object)
  416.    OBJECT *Object;
  417. {
  418.    QUARTIC *Shape = (QUARTIC *) Object;
  419.    int i;
  420.    for (i=0;i<35;i++)
  421.       Shape->Coeffs[i] *= -1.0;
  422. }
  423.  
  424. /* Given a set of power values from a term in the general quartic, return
  425.    the index of the term. */
  426. int
  427. roll_term_indices(i, j, k)
  428.    int i, j, k;
  429. {
  430.    int l, t;
  431.    l = 4 - (i + j + k);
  432.    if (i < 0 || i > 4 || j < 0 || j > 4 ||
  433.        k < 0 || k > 4 || l < 0 || l > 4) {
  434.       printf("Error in 'roll_term_indices': i = %d, j = %d, k = %d, l = %d\n",
  435.          i, j, k, l);
  436.       exit(1);
  437.       }
  438.    for (t=0;t<term_counts[4];t++)
  439.       if (i == terms4[t][0] && j == terms4[t][1] &&
  440.           k == terms4[t][2] && l == terms4[t][3])
  441.      return t;
  442.    printf("Term not  in 'roll_term_indices': i = %d, j = %d, k = %d, l = %d\n",
  443.          i, j, k, l);
  444.    exit(1);
  445. }
  446.  
  447. /* Given the total power and the index, return the individual powers */
  448. void
  449. unroll_term_indices( Power, index, i, j, k, l)
  450.    int Power, index, *i, *j, *k, *l;
  451. {
  452.    if (Power == 4) {
  453.       *i = terms4[index][0];
  454.       *j = terms4[index][1];
  455.       *k = terms4[index][2];
  456.       *l = terms4[index][3];
  457.       }
  458.    else if (Power == 3) {
  459.       *i = terms3[index][0];
  460.       *j = terms3[index][1];
  461.       *k = terms3[index][2];
  462.       *l = terms3[index][3];
  463.       }
  464.    else if (Power == 2) {
  465.       *i = terms2[index][0];
  466.       *j = terms2[index][1];
  467.       *k = terms2[index][2];
  468.       *l = terms2[index][3];
  469.       }
  470.    else if (Power == 1) {
  471.       *i = terms1[index][0];
  472.       *j = terms1[index][1];
  473.       *k = terms1[index][2];
  474.       *l = terms1[index][3];
  475.       }
  476.    else {
  477.       *i = 0;
  478.       *j = 0;
  479.       *k = 0;
  480.       *l = 0;
  481.       }
  482. }
  483.  
  484. DBL
  485. do_partial_term(q, row, Power, i, j, k, l)
  486.    MATRIX *q;
  487.    int row, Power, i, j, k, l;
  488. {
  489.    DBL result;
  490.    result = ((DBL) factorials[Power] /
  491.           (factorials[i]*factorials[j]*factorials[k]*factorials[l]));
  492.    if (i>0) result *= POW((*q)[0][row], (DBL) i);
  493.    if (j>0) result *= POW((*q)[1][row], (DBL) j);
  494.    if (k>0) result *= POW((*q)[2][row], (DBL) k);
  495.    if (l>0) result *= POW((*q)[3][row], (DBL) l);
  496.    return result;
  497. }
  498.  
  499. /* Using the transformation matrix q, transform the general quartic equation
  500.    given by a. */
  501. void
  502. Transform_Quartic(Shape, q)
  503.    QUARTIC *Shape;
  504.    MATRIX *q;
  505. {
  506.    int term_index, partial_index;
  507.    int ip, i, i0, i1, i2, i3;
  508.    int jp, j, j0, j1, j2, j3;
  509.    int kp, k, k0, k1, k2, k3;
  510.    DBL b[35], partial_term;
  511.    DBL tempx, tempy, tempz;
  512.  
  513.    /* Trim out any really low values from the transform. */
  514.    for (i=0;i<4;i++) for (j=0;j<4;j++)
  515.       if ((*q)[i][j] > -EPSILON && (*q)[i][j] < EPSILON) (*q)[i][j] = 0.0;
  516.  
  517.    for (i=0;i<35;i++) b[i] = 0.0;
  518.  
  519.    for (term_index=0;term_index<term_counts[4];term_index++) {
  520.       if (Shape->Coeffs[term_index] == 0.0) continue;
  521.       ip = terms4[term_index][0];    /* Power of original x term */
  522.       jp = terms4[term_index][1];    /* Power of original y term */
  523.       kp = terms4[term_index][2];    /* Power of original z term */
  524.       /* Step through terms in: (q[0][0]*x+q[0][1]*y+q[0][2]*z+q[0][3])^i */
  525.       for (i=0;i<term_counts[ip];i++) {
  526.      unroll_term_indices(ip, i, &i0, &i1, &i2, &i3);
  527.      tempx = do_partial_term(q, 0, ip, i0, i1, i2, i3);
  528.      if (tempx == 0.0) continue;
  529.  
  530.       /* Step through terms in: (q[1][0]*x+q[1][1]*y+q[1][2]*z+q[1][3])^j */
  531.      for (j=0;j<term_counts[jp];j++) {
  532.         unroll_term_indices(jp, j, &j0, &j1, &j2, &j3);
  533.         tempy = do_partial_term(q, 1, jp, j0, j1, j2, j3);
  534.         if (tempy == 0.0) continue;
  535.  
  536.       /* Step through terms in: (q[2][0]*x+q[2][1]*y+q[2][2]*z+q[2][3])^k */
  537.         for (k=0;k<term_counts[kp];k++) {
  538.            unroll_term_indices(kp, k, &k0, &k1, &k2, &k3);
  539.            tempz = do_partial_term(q, 2, kp, k0, k1, k2, k3);
  540.             if (tempz == 0.0) continue;
  541.  
  542.            /* Figure out it's index, and add it into the result */
  543.            partial_index = roll_term_indices(i0+j0+k0,i1+j1+k1,i2+j2+k2);
  544.            partial_term = Shape->Coeffs[term_index] * tempx * tempy * tempz;
  545.            b[partial_index] += partial_term;
  546.            }
  547.         }
  548.      }
  549.       }
  550.    for (i=0;i<35;i++) {
  551.       if (b[i] > -EPSILON && b[i] < EPSILON)
  552.      Shape->Coeffs[i] = 0.0;
  553.       else
  554.      Shape->Coeffs[i] = b[i];
  555.       }
  556. }
  557.