home *** CD-ROM | disk | FTP | other *** search
/ Chestnut's Multimedia Mania / MM_MANIA.ISO / graphics / povsrc20 / quadrics.c < prev    next >
C/C++ Source or Header  |  1993-07-29  |  10KB  |  356 lines
  1. /****************************************************************************
  2. *                quadrics.c
  3. *
  4. *  This module implements the code for the quadric shape primitive.
  5. *
  6. *  from Persistence of Vision Raytracer
  7. *  Copyright 1993 Persistence of Vision Team
  8. *---------------------------------------------------------------------------
  9. *  NOTICE: This source code file is provided so that users may experiment
  10. *  with enhancements to POV-Ray and to port the software to platforms other 
  11. *  than those supported by the POV-Ray Team.  There are strict rules under
  12. *  which you are permitted to use this file.  The rules are in the file
  13. *  named POVLEGAL.DOC which should be distributed with this file. If 
  14. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  15. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  16. *  Forum.  The latest version of POV-Ray may be found there as well.
  17. *
  18. * This program is based on the popular DKB raytracer version 2.12.
  19. * DKBTrace was originally written by David K. Buck.
  20. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  21. *
  22. *****************************************************************************/
  23.  
  24. #include "frame.h"
  25. #include "vector.h"
  26. #include "povproto.h"
  27.  
  28. METHODS Quadric_Methods =
  29.   { 
  30.   All_Quadric_Intersections,
  31.   Inside_Quadric, Quadric_Normal,
  32.   Copy_Quadric,
  33.   Translate_Quadric, Rotate_Quadric,
  34.   Scale_Quadric, Transform_Quadric, Invert_Quadric,
  35.   Destroy_Quadric
  36. };
  37.  
  38. extern RAY *CM_Ray;
  39. extern long Ray_Quadric_Tests, Ray_Quadric_Tests_Succeeded;
  40.  
  41. int All_Quadric_Intersections (Object, Ray, Depth_Stack)
  42. OBJECT *Object;
  43. RAY *Ray;
  44. ISTACK *Depth_Stack;
  45.   {
  46.   DBL Depth1, Depth2;
  47.   VECTOR IPoint;
  48.   register int Intersection_Found;
  49.  
  50.   Intersection_Found = FALSE;
  51.  
  52.   if (Intersect_Quadric (Ray, (QUADRIC *) Object, &Depth1, &Depth2))
  53.     {
  54.     VScale (IPoint, Ray->Direction, Depth1);
  55.     VAddEq (IPoint, Ray->Initial);
  56.  
  57.     if (Point_In_Clip (&IPoint, Object->Clip))
  58.       {
  59.       push_entry(Depth1,IPoint,Object,Depth_Stack);
  60.       Intersection_Found = TRUE;
  61.       }
  62.  
  63.     if (Depth2 != Depth1)
  64.       {
  65.       VScale (IPoint, Ray->Direction, Depth2);
  66.       VAddEq (IPoint, Ray->Initial);
  67.  
  68.       if (Point_In_Clip (&IPoint, Object->Clip))
  69.         {
  70.         push_entry(Depth2,IPoint,Object,Depth_Stack);
  71.         Intersection_Found = TRUE;
  72.         }
  73.       }
  74.     }
  75.   return (Intersection_Found);
  76.   }
  77.  
  78. int Intersect_Quadric (Ray, Quadric, Depth1, Depth2)
  79. RAY *Ray;
  80. QUADRIC *Quadric;
  81. DBL *Depth1, *Depth2;
  82.   {
  83.   register DBL Square_Term, Linear_Term, Constant_Term, Temp_Term;
  84.   register DBL Determinant, Determinant_2, A2, BMinus;
  85.  
  86.   Ray_Quadric_Tests++;
  87.   if (!Ray->Quadric_Constants_Cached)
  88.     Make_Ray(Ray);
  89.  
  90.   if (Quadric->Non_Zero_Square_Term)
  91.     {
  92.     VDot (Square_Term, Quadric->Square_Terms, Ray->Direction_2);
  93.     VDot (Temp_Term, Quadric->Mixed_Terms, Ray->Mixed_Dir_Dir);
  94.     Square_Term += Temp_Term;
  95.     }
  96.   else
  97.     Square_Term = 0.0;
  98.  
  99.   VDot (Linear_Term, Quadric->Square_Terms, Ray->Initial_Direction);
  100.   Linear_Term *= 2.0;
  101.   VDot (Temp_Term, Quadric->Terms, Ray->Direction);
  102.   Linear_Term += Temp_Term;
  103.   VDot (Temp_Term, Quadric->Mixed_Terms, Ray->Mixed_Init_Dir);
  104.   Linear_Term += Temp_Term;
  105.  
  106.   if (Ray == CM_Ray)
  107.     if (!Quadric->Constant_Cached)
  108.       {
  109.       VDot (Constant_Term, Quadric->Square_Terms, Ray->Initial_2);
  110.       VDot (Temp_Term, Quadric->Terms, Ray->Initial);
  111.       Constant_Term +=  Temp_Term + Quadric->Constant;
  112.       Quadric->CM_Constant = Constant_Term;
  113.       Quadric->Constant_Cached = TRUE;
  114.       }
  115.     else
  116.       Constant_Term = Quadric->CM_Constant;
  117.   else
  118.     {
  119.     VDot (Constant_Term, Quadric->Square_Terms, Ray->Initial_2);
  120.     VDot (Temp_Term, Quadric->Terms, Ray->Initial);
  121.     Constant_Term += Temp_Term + Quadric->Constant;
  122.     }
  123.  
  124.   VDot (Temp_Term, Quadric->Mixed_Terms, 
  125.     Ray->Mixed_Initial_Initial);
  126.   Constant_Term += Temp_Term;
  127.  
  128.   if (Square_Term != 0.0)
  129.     {
  130.     /* The equation is quadratic - find its roots */
  131.  
  132.     Determinant_2 = Linear_Term * Linear_Term - 4.0 * Square_Term * Constant_Term;
  133.  
  134.     if (Determinant_2 < 0.0)
  135.       return (FALSE);
  136.  
  137.     Determinant = sqrt (Determinant_2);
  138.     A2 = Square_Term * 2.0;
  139.     BMinus = Linear_Term * -1.0;
  140.  
  141.     *Depth1 = (BMinus + Determinant) / A2;
  142.     *Depth2 = (BMinus - Determinant) / A2;
  143.     }
  144.   else
  145.     {
  146.     /* There are no quadratic terms.  Solve the linear equation instead. */
  147.     if (Linear_Term == 0.0)
  148.       return (FALSE);
  149.  
  150.     *Depth1 = Constant_Term * -1.0 / Linear_Term;
  151.     *Depth2 = *Depth1;
  152.     }
  153.  
  154.   if ((*Depth1 < Small_Tolerance) || (*Depth1 > Max_Distance))
  155.     if ((*Depth2 < Small_Tolerance) || (*Depth2 > Max_Distance))
  156.       return (FALSE);
  157.     else
  158.       *Depth1 = *Depth2;
  159.   else
  160.     if ((*Depth2 < Small_Tolerance) || (*Depth2 > Max_Distance))
  161.       *Depth2 = *Depth1;
  162.  
  163.   Ray_Quadric_Tests_Succeeded++;
  164.   return (TRUE);
  165.   }
  166.  
  167. int Inside_Quadric (IPoint, Object)
  168. VECTOR *IPoint;
  169. OBJECT *Object;
  170.   {
  171.   QUADRIC *Quadric = (QUADRIC *) Object;
  172.   VECTOR New_Point;
  173.   register DBL Result, Linear_Term, Square_Term;
  174.  
  175.   VDot (Linear_Term, *IPoint, Quadric->Terms);
  176.   Result = Linear_Term + Quadric->Constant;
  177.   VSquareTerms (New_Point, *IPoint);
  178.   VDot (Square_Term, New_Point, Quadric->Square_Terms);
  179.   Result += Square_Term;
  180.   Result += Quadric->Mixed_Terms.x * (IPoint->x) * (IPoint->y)
  181.     + Quadric->Mixed_Terms.y * (IPoint->x) * (IPoint->z)
  182.       + Quadric->Mixed_Terms.z * (IPoint->y) * (IPoint->z);
  183.  
  184.   if (Result < Small_Tolerance)
  185.     return (TRUE);
  186.  
  187.   return (FALSE);
  188.   }
  189.  
  190. void Quadric_Normal (Result, Object, IPoint)
  191. VECTOR *Result, *IPoint;
  192. OBJECT *Object;
  193.   {
  194.   QUADRIC *Intersection_Quadric = (QUADRIC *) Object;
  195.   VECTOR Derivative_Linear;
  196.   DBL Len;
  197.  
  198.   VScale (Derivative_Linear, Intersection_Quadric->Square_Terms, 2.0);
  199.   VEvaluate (*Result, Derivative_Linear, *IPoint);
  200.   VAdd (*Result, *Result, Intersection_Quadric->Terms);
  201.  
  202.   Result->x += 
  203.   Intersection_Quadric->Mixed_Terms.x * IPoint->y +
  204.   Intersection_Quadric->Mixed_Terms.y * IPoint->z;
  205.  
  206.  
  207.   Result->y +=
  208.   Intersection_Quadric->Mixed_Terms.x * IPoint->x +
  209.   Intersection_Quadric->Mixed_Terms.z * IPoint->z;
  210.  
  211.   Result->z +=
  212.   Intersection_Quadric->Mixed_Terms.y * IPoint->x +
  213.   Intersection_Quadric->Mixed_Terms.z * IPoint->y;
  214.  
  215.   Len = Result->x * Result->x + Result->y * Result->y + Result->z * Result->z;
  216.   Len = sqrt(Len);
  217.   if (Len == 0.0) 
  218.     {
  219.     /* The normal is not defined at this point of the surface.  Set it
  220.          to any arbitrary direction. */
  221.     Result->x = 1.0;
  222.     Result->y = 0.0;
  223.     Result->z = 0.0;
  224.     }
  225.   else 
  226.     {
  227.     Result->x /= Len;        /* normalize the normal */
  228.     Result->y /= Len;
  229.     Result->z /= Len;
  230.     }
  231.   }
  232.  
  233.   void Transform_Quadric (Object, Trans)
  234.     OBJECT *Object;
  235. TRANSFORM *Trans;
  236.   {
  237.   QUADRIC *Quadric=(QUADRIC *)Object;
  238.   MATRIX Quadric_Matrix, Transform_Transposed;
  239.  
  240.   Quadric_To_Matrix (Quadric, (MATRIX *) &Quadric_Matrix[0][0]);
  241.   MTimes ((MATRIX *) &Quadric_Matrix[0][0], (MATRIX *) &(Trans->inverse[0][0]), (MATRIX *) &Quadric_Matrix[0][0]);
  242.   MTranspose ((MATRIX *) &Transform_Transposed[0][0], (MATRIX *) &(Trans->inverse[0][0]));
  243.   MTimes ((MATRIX *) &Quadric_Matrix[0][0], (MATRIX *) &Quadric_Matrix[0][0], (MATRIX *) &Transform_Transposed[0][0]);
  244.   Matrix_To_Quadric ((MATRIX *) &Quadric_Matrix[0][0], Quadric);
  245.   }
  246.  
  247. void Quadric_To_Matrix (Quadric, Matrix)
  248. QUADRIC *Quadric;
  249. MATRIX *Matrix;
  250.   {
  251.   MZero (Matrix);
  252.   (*Matrix)[0][0] = Quadric->Square_Terms.x;
  253.   (*Matrix)[1][1] = Quadric->Square_Terms.y;
  254.   (*Matrix)[2][2] = Quadric->Square_Terms.z;
  255.   (*Matrix)[0][1] = Quadric->Mixed_Terms.x;
  256.   (*Matrix)[0][2] = Quadric->Mixed_Terms.y;
  257.   (*Matrix)[0][3] = Quadric->Terms.x;
  258.   (*Matrix)[1][2] = Quadric->Mixed_Terms.z;
  259.   (*Matrix)[1][3] = Quadric->Terms.y;
  260.   (*Matrix)[2][3] = Quadric->Terms.z;
  261.   (*Matrix)[3][3] = Quadric->Constant;
  262.   }
  263.  
  264. void Matrix_To_Quadric (Matrix, Quadric)
  265. MATRIX *Matrix;
  266. QUADRIC *Quadric;
  267.   {
  268.   Quadric->Square_Terms.x = (*Matrix)[0][0];
  269.   Quadric->Square_Terms.y = (*Matrix)[1][1];
  270.   Quadric->Square_Terms.z = (*Matrix)[2][2];
  271.   Quadric->Mixed_Terms.x = (*Matrix)[0][1] + (*Matrix)[1][0];
  272.   Quadric->Mixed_Terms.y = (*Matrix)[0][2] + (*Matrix)[2][0];
  273.   Quadric->Terms.x = (*Matrix)[0][3] + (*Matrix)[3][0];
  274.   Quadric->Mixed_Terms.z = (*Matrix)[1][2] + (*Matrix)[2][1];
  275.   Quadric->Terms.y = (*Matrix)[1][3] + (*Matrix)[3][1];
  276.   Quadric->Terms.z = (*Matrix)[2][3] + (*Matrix)[3][2];
  277.   Quadric->Constant = (*Matrix)[3][3];
  278.   }
  279.  
  280. void Translate_Quadric (Object, Vector)
  281. OBJECT *Object;
  282. VECTOR *Vector;
  283.   {
  284.   TRANSFORM Trans;
  285.  
  286.   Compute_Translation_Transform (&Trans, Vector);
  287.   Transform_Quadric (Object, &Trans);
  288.  
  289.   }
  290.  
  291. void Rotate_Quadric (Object, Vector)
  292. OBJECT *Object;
  293. VECTOR *Vector;
  294.   {
  295.   TRANSFORM Trans;
  296.  
  297.   Compute_Rotation_Transform (&Trans, Vector);
  298.   Transform_Quadric (Object, &Trans);
  299.   }
  300.  
  301. void Scale_Quadric (Object, Vector)
  302. OBJECT *Object;
  303. VECTOR *Vector;
  304.   {
  305.   TRANSFORM Trans;
  306.  
  307.   Compute_Scaling_Transform (&Trans, Vector);
  308.   Transform_Quadric (Object, &Trans);
  309.   }
  310.  
  311. void Invert_Quadric (Object)
  312. OBJECT *Object;
  313.   {
  314.   QUADRIC *Quadric = (QUADRIC *) Object;
  315.  
  316.   VScaleEq (Quadric->Square_Terms, -1.0);
  317.   VScaleEq (Quadric->Mixed_Terms, -1.0);
  318.   VScaleEq (Quadric->Terms, -1.0);
  319.   Quadric->Constant *= -1.0;
  320.   }
  321.  
  322. QUADRIC *Create_Quadric()
  323.   {
  324.   QUADRIC *New;
  325.  
  326.   if ((New = (QUADRIC *) malloc (sizeof (QUADRIC))) == NULL)
  327.     MAError ("quadric");
  328.  
  329.   INIT_OBJECT_FIELDS(New, QUADRIC_OBJECT, &Quadric_Methods)
  330.     Make_Vector (&(New->Square_Terms), 1.0, 1.0, 1.0);
  331.   Make_Vector (&(New->Mixed_Terms), 0.0, 0.0, 0.0);
  332.   Make_Vector (&(New->Terms), 0.0, 0.0, 0.0);
  333.   New->Constant = 1.0;
  334.   New->CM_Constant = HUGE_VAL;
  335.   New->Constant_Cached = FALSE;
  336.   New->Non_Zero_Square_Term = FALSE;
  337.   return (New);
  338.   }
  339.  
  340. void *Copy_Quadric (Object)
  341. OBJECT *Object;
  342.   {
  343.   QUADRIC *New;
  344.  
  345.   New = Create_Quadric ();
  346.   *New = *((QUADRIC *) Object);
  347.  
  348.   return (New);
  349.   }
  350.  
  351. void Destroy_Quadric (Object)
  352. OBJECT *Object;
  353.   {
  354.   free (Object);
  355.   }
  356.