/*****************************************************************************
*
*                                   quadrics.c
*
*   from DKBTrace (c) 1990  David Buck
*
*  This module implements the code for the quadric shape primitive.
*
* This software is freely distributable. The source and/or object code may be
* copied or uploaded to communications services so long as this notice remains
* at the top of each file.  If any changes are made to the program, you must
* clearly indicate in the documentation and in the programs startup message
* who it was who made the changes. The documentation should also describe what
* those changes were. This software may not be included in whole or in
* part into any commercial package without the express written consent of the
* author.  It may, however, be included in other public domain or freely
* distributed software so long as the proper credit for the software is given.
*
* This software is provided as is without any guarantees or warranty. Although
* the author has attempted to find and correct any bugs in the software, he
* is not responsible for any damage caused by the use of the software.  The
* author is under no obligation to provide service, corrections, or upgrades
* to this package.
*
* Despite all the legal stuff above, if you do find bugs, I would like to hear
* about them.  Also, if you have any comments or questions, you may contact me
* at the following address:
*
*     David Buck
*     22C Sonnet Cres.
*     Nepean Ontario
*     Canada, K2H 8W7
*
*  I can also be reached on the following bulleton boards:
*
*     OMX              (613) 731-3419
*     Mystic           (613) 596-4249  or  (613) 596-4772
*
*  Fidonet:   1:163/109.9
*  Internet:  dbuck@ccs.carleton.ca
*  The "You Can Call Me RAY" BBS    (708) 358-5611
*
*  IBM Port by Aaron A. Collins. Aaron may be reached on the following BBS'es:
*
*     The "You Can Call Me RAY" BBS (708) 358-5611
*     The Information Exchange BBS  (708) 945-5575
*
*****************************************************************************/


#include "frame.h"
#include "vector.h"
#include "dkbproto.h"

METHODS Quadric_Methods =
   { Object_Intersect, All_Quadric_Intersections,
     Inside_Quadric, Quadric_Normal,
     Copy_Quadric,
     Translate_Quadric, Rotate_Quadric,
     Scale_Quadric, Invert_Quadric};

extern RAY *VP_Ray;
extern long Ray_Quadric_Tests, Ray_Quadric_Tests_Succeeded;

int All_Quadric_Intersections (Object, Ray, Depth_Queue)
   OBJECT *Object;
   RAY *Ray;
   PRIOQ *Depth_Queue;
   {
   QUADRIC *Shape = (QUADRIC *) Object;
   DBL Depth1, Depth2;
   VECTOR Intersection_Point;
   INTERSECTION Local_Element;
   register int Intersection_Found;

   Intersection_Found = FALSE;
   if (Intersect_Quadric (Ray, Shape, &Depth1, &Depth2))
      {
      Local_Element.Depth = Depth1;
      Local_Element.Object = Shape -> Parent_Object;
      VScale (Intersection_Point, Ray -> Direction, Depth1);
      VAdd (Intersection_Point, Intersection_Point, Ray -> Initial);
      Local_Element.Point = Intersection_Point;
      Local_Element.Shape = (SHAPE *)Shape;
      pq_add (Depth_Queue, &Local_Element);
      Intersection_Found = TRUE;

      if (Depth2 != Depth1)
         {
         Local_Element.Depth = Depth2;
         Local_Element.Object = Shape -> Parent_Object;
         VScale (Intersection_Point, Ray -> Direction, Depth2);
         VAdd (Intersection_Point, Intersection_Point, Ray -> Initial);
         Local_Element.Point = Intersection_Point;
         Local_Element.Shape = (SHAPE *) Shape;
         pq_add (Depth_Queue, &Local_Element);
         Intersection_Found = TRUE;
         }
      }
   return (Intersection_Found);
   }

int Intersect_Quadric (Ray, Shape, Depth1, Depth2)
   RAY *Ray;
   QUADRIC *Shape;
   DBL *Depth1, *Depth2;
   {
   register DBL Square_Term, Linear_Term, Constant_Term, Temp_Term;
   register DBL Determinant, Determinant_2, A2, BMinus;

   Ray_Quadric_Tests++;
   if (!Ray->Quadric_Constants_Cached)
      Make_Ray(Ray);

   if (Shape -> Non_Zero_Square_Term)
      {
      VDot (Square_Term, Shape -> Object_2_Terms, Ray -> Direction_2);
      VDot (Temp_Term, Shape -> Object_Mixed_Terms, Ray -> Mixed_Dir_Dir);
      Square_Term += Temp_Term;
      }
   else
      Square_Term = 0.0;

   VDot (Linear_Term, Shape -> Object_2_Terms, Ray -> Initial_Direction);
   Linear_Term *= 2.0;
   VDot (Temp_Term, Shape -> Object_Terms, Ray -> Direction);
   Linear_Term += Temp_Term;
   VDot (Temp_Term, Shape -> Object_Mixed_Terms, Ray -> Mixed_Init_Dir);
   Linear_Term += Temp_Term;

   if (Ray == VP_Ray)
      if (!Shape -> Constant_Cached)
         {
         VDot (Constant_Term, Shape -> Object_2_Terms, Ray -> Initial_2);
         VDot (Temp_Term, Shape -> Object_Terms, Ray -> Initial);
         Constant_Term +=  Temp_Term + Shape -> Object_Constant;
         Shape -> Object_VP_Constant = Constant_Term;
         Shape -> Constant_Cached = TRUE;
         }
      else
         Constant_Term = Shape -> Object_VP_Constant;
   else
      {
      VDot (Constant_Term, Shape -> Object_2_Terms, Ray -> Initial_2);
      VDot (Temp_Term, Shape -> Object_Terms, Ray -> Initial);
      Constant_Term += Temp_Term + Shape -> Object_Constant;
      }

   VDot (Temp_Term, Shape -> Object_Mixed_Terms, 
         Ray -> Mixed_Initial_Initial);
   Constant_Term += Temp_Term;

   if (Square_Term != 0.0)
      {
  /* The equation is quadratic - find its roots */

      Determinant_2 = Linear_Term * Linear_Term - 4.0 * Square_Term * Constant_Term;

      if (Determinant_2 < 0.0)
         return (FALSE);

      Determinant = sqrt (Determinant_2);
      A2 = Square_Term * 2.0;
      BMinus = Linear_Term * -1.0;

      *Depth1 = (BMinus + Determinant) / A2;
      *Depth2 = (BMinus - Determinant) / A2;
      }
   else
      {
    /* There are no quadratic terms.  Solve the linear equation instead. */
      if (Linear_Term == 0.0)
         return (FALSE);

      *Depth1 = Constant_Term * -1.0 / Linear_Term;
      *Depth2 = *Depth1;
      }

   if ((*Depth1 < Small_Tolerance) || (*Depth1 > Max_Distance))
      if ((*Depth2 < Small_Tolerance) || (*Depth2 > Max_Distance))
         return (FALSE);
      else
         *Depth1 = *Depth2;
   else
      if ((*Depth2 < Small_Tolerance) || (*Depth2 > Max_Distance))
        *Depth2 = *Depth1;

   Ray_Quadric_Tests_Succeeded++;
   return (TRUE);
   }

int Inside_Quadric (Test_Point, Object)
   VECTOR *Test_Point;
   OBJECT *Object;
   {
   QUADRIC *Shape = (QUADRIC *) Object;
   VECTOR New_Point;
   register DBL Result, Linear_Term, Square_Term;

   VDot (Linear_Term, *Test_Point, Shape -> Object_Terms);
   Result = Linear_Term + Shape -> Object_Constant;
   VSquareTerms (New_Point, *Test_Point);
   VDot (Square_Term, New_Point, Shape -> Object_2_Terms);
   Result += Square_Term;
   Result += Shape -> Object_Mixed_Terms.x * (Test_Point -> x) * (Test_Point -> y)
            + Shape -> Object_Mixed_Terms.y * (Test_Point -> x) * (Test_Point -> z)
            + Shape -> Object_Mixed_Terms.z * (Test_Point -> y) * (Test_Point -> z);

   if (Result < Small_Tolerance)
      return (TRUE);

   return (FALSE);
   }

void Quadric_Normal (Result, Object, Intersection_Point)
   VECTOR *Result, *Intersection_Point;
   OBJECT *Object;
   {
   QUADRIC *Intersection_Shape = (QUADRIC *) Object;
   VECTOR Derivative_Linear;
   DBL Len;

   VScale (Derivative_Linear, Intersection_Shape -> Object_2_Terms, 2.0);
   VEvaluate (*Result, Derivative_Linear, *Intersection_Point);
   VAdd (*Result, *Result, Intersection_Shape -> Object_Terms);

   Result -> x += 
         Intersection_Shape->Object_Mixed_Terms.x * Intersection_Point->y +
         Intersection_Shape->Object_Mixed_Terms.y * Intersection_Point->z;


   Result -> y +=
         Intersection_Shape->Object_Mixed_Terms.x * Intersection_Point->x +
         Intersection_Shape->Object_Mixed_Terms.z * Intersection_Point->z;

   Result -> z +=
         Intersection_Shape->Object_Mixed_Terms.y * Intersection_Point->x +
         Intersection_Shape->Object_Mixed_Terms.z * Intersection_Point->y;

   Len = Result->x * Result->x + Result->y * Result->y + Result->z * Result->z;
   Len = sqrt(Len);
   if (Len == 0.0) {
      /* The normal is not defined at this point of the surface.  Set it
         to any arbitrary direction. */
      Result->x = 1.0;
      Result->y = 0.0;
      Result->z = 0.0;
      }
   else {
      Result->x /= Len;		/* normalize the normal */
      Result->y /= Len;
      Result->z /= Len;
      }
   }

void *Copy_Quadric (Object)
   OBJECT *Object;
   {
   QUADRIC *New_Shape;

   New_Shape = Get_Quadric_Shape ();
   *New_Shape = *((QUADRIC *) Object);
   New_Shape -> Next_Object = NULL;

   if (New_Shape->Shape_Texture != NULL)
      New_Shape->Shape_Texture = Copy_Texture (New_Shape->Shape_Texture);

   return (New_Shape);
   }

void Transform_Quadric (Shape, Transformation)
   QUADRIC *Shape;
   TRANSFORMATION *Transformation;
   {
   MATRIX Quadric_Matrix, Transform_Transposed;

   Quadric_To_Matrix (Shape, (MATRIX *) &Quadric_Matrix[0][0]);
   MTimes ((MATRIX *) &Quadric_Matrix[0][0], (MATRIX *) &(Transformation -> inverse[0][0]), (MATRIX *) &Quadric_Matrix[0][0]);
   MTranspose ((MATRIX *) &Transform_Transposed[0][0], (MATRIX *) &(Transformation -> inverse[0][0]));
   MTimes ((MATRIX *) &Quadric_Matrix[0][0], (MATRIX *) &Quadric_Matrix[0][0], (MATRIX *) &Transform_Transposed[0][0]);
   Matrix_To_Quadric ((MATRIX *) &Quadric_Matrix[0][0], Shape);
   }

void Quadric_To_Matrix (Quadric, Matrix)
   QUADRIC *Quadric;
   MATRIX *Matrix;
   {
   MZero (Matrix);
   (*Matrix)[0][0] = Quadric -> Object_2_Terms.x;
   (*Matrix)[1][1] = Quadric -> Object_2_Terms.y;
   (*Matrix)[2][2] = Quadric -> Object_2_Terms.z;
   (*Matrix)[0][1] = Quadric -> Object_Mixed_Terms.x;
   (*Matrix)[0][2] = Quadric -> Object_Mixed_Terms.y;
   (*Matrix)[0][3] = Quadric -> Object_Terms.x;
   (*Matrix)[1][2] = Quadric -> Object_Mixed_Terms.z;
   (*Matrix)[1][3] = Quadric -> Object_Terms.y;
   (*Matrix)[2][3] = Quadric -> Object_Terms.z;
   (*Matrix)[3][3] = Quadric -> Object_Constant;
   }

void Matrix_To_Quadric (Matrix, Quadric)
   MATRIX *Matrix;
   QUADRIC *Quadric;
   {
   Quadric -> Object_2_Terms.x = (*Matrix)[0][0];
   Quadric -> Object_2_Terms.y = (*Matrix)[1][1];
   Quadric -> Object_2_Terms.z = (*Matrix)[2][2];
   Quadric -> Object_Mixed_Terms.x = (*Matrix)[0][1] + (*Matrix)[1][0];
   Quadric -> Object_Mixed_Terms.y = (*Matrix)[0][2] + (*Matrix)[2][0];
   Quadric -> Object_Terms.x = (*Matrix)[0][3] + (*Matrix)[3][0];
   Quadric -> Object_Mixed_Terms.z = (*Matrix)[1][2] + (*Matrix)[2][1];
   Quadric -> Object_Terms.y = (*Matrix)[1][3] + (*Matrix)[3][1];
   Quadric -> Object_Terms.z = (*Matrix)[2][3] + (*Matrix)[3][2];
   Quadric -> Object_Constant = (*Matrix)[3][3];
   }

void Translate_Quadric (Object, Vector)
   OBJECT *Object;
   VECTOR *Vector;
   {
   TRANSFORMATION Transformation;

   Get_Translation_Transformation (&Transformation, Vector);
   Transform_Quadric ((QUADRIC *) Object, &Transformation);

   Translate_Texture (&((QUADRIC *) Object)->Shape_Texture, Vector);
   }

void Rotate_Quadric (Object, Vector)
   OBJECT *Object;
   VECTOR *Vector;
   {
   TRANSFORMATION Transformation;

   Get_Rotation_Transformation (&Transformation, Vector);
   Transform_Quadric ((QUADRIC *) Object, &Transformation);

   Rotate_Texture (&((QUADRIC *) Object)->Shape_Texture, Vector);
   }

void Scale_Quadric (Object, Vector)
   OBJECT *Object;
   VECTOR *Vector;
   {
   TRANSFORMATION Transformation;

   Get_Scaling_Transformation (&Transformation, Vector);
   Transform_Quadric ((QUADRIC *) Object, &Transformation);

   Scale_Texture (&((QUADRIC *) Object)->Shape_Texture, Vector);
   }

void Invert_Quadric (Object)
   OBJECT *Object;
   {
   QUADRIC *Shape = (QUADRIC *) Object;

   VScale (Shape -> Object_2_Terms, Shape -> Object_2_Terms, -1.0);
   VScale (Shape -> Object_Mixed_Terms, Shape -> Object_Mixed_Terms, -1.0);
   VScale (Shape -> Object_Terms, Shape -> Object_Terms, -1.0);
   Shape -> Object_Constant *= -1.0;
   }