/****************************************************************************
*
*  ATTENTION!!!
*
*  THIS FILE HAS BEEN MODIFIED!!! IT IS NOT PART OF THE OFFICAL
*  POV-RAY 2.2 DISTRIBUTION!!!
*
*  THIS FILE IS PART OF "FASTER THAN POV-RAY" (VERSION 1.1),
*  A SPED-UP VERSION OF POV-RAY 2.2. USE AT YOUR OWN RISK!!!!!!
*
*  New files: addon0.c, addon1.c, addon2.c, addon3.c, addon.h
*
*  The additional modules were written by Dieter Bayer.
*
*  Send comments, suggestions, bugs, ideas ... to:
*
*  dieter@cip.e-technik.uni-erlangen.de
*
*  All changed/added lines are enclosed in #ifdef DB_CODE ... #endif
*
*  The vista projection was taken from:
*
*    A. Hashimoto, T. Akimoto, K. Mase, and Y. Suenaga, 
*    "Vista Ray-Tracing: High Speed Ray Tracing Using Perspective
*    Projection Image", New Advances in Computer Graphics, Proceedings
*    of CG International '89, R. A. Earnshaw, B. Wyvill (Eds.), 
*    Springer, ..., pp. 549-560
*
*  The idea for the light buffer was taken from:
*
*    E. Haines and D. Greenberg, "The Light Buffer: A Shadow-Testing 
*    Accelerator", IEEE CG&A, Vol. 6, No. 9, Sept. 1986, pp. 6-16
*
*****************************************************************************/

/****************************************************************************
*                quadrics.c
*
*  This module implements the code for the quadric shape primitive.
*
*  from Persistence of Vision Raytracer
*  Copyright 1993 Persistence of Vision Team
*---------------------------------------------------------------------------
*  NOTICE: This source code file is provided so that users may experiment
*  with enhancements to POV-Ray and to port the software to platforms other 
*  than those supported by the POV-Ray Team.  There are strict rules under
*  which you are permitted to use this file.  The rules are in the file
*  named POVLEGAL.DOC which should be distributed with this file. If 
*  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
*  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
*  Forum.  The latest version of POV-Ray may be found there as well.
*
* This program is based on the popular DKB raytracer version 2.12.
* DKBTrace was originally written by David K. Buck.
* DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
*
*****************************************************************************/

#include "frame.h"
#include "vector.h"
#include "povproto.h"

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

extern RAY *CM_Ray;
extern long Ray_Quadric_Tests, Ray_Quadric_Tests_Succeeded;

int All_Quadric_Intersections (Object, Ray, Depth_Stack)
OBJECT *Object;
RAY *Ray;
ISTACK *Depth_Stack;
  {
  DBL Depth1, Depth2;
  VECTOR IPoint;
  register int Intersection_Found;

  Intersection_Found = FALSE;

  if (Intersect_Quadric (Ray, (QUADRIC *) Object, &Depth1, &Depth2))
    {
    VScale (IPoint, Ray->Direction, Depth1);
    VAddEq (IPoint, Ray->Initial);

    if (Point_In_Clip (&IPoint, Object->Clip))
      {
      push_entry(Depth1,IPoint,Object,Depth_Stack);
      Intersection_Found = TRUE;
      }

    if (Depth2 != Depth1)
      {
      VScale (IPoint, Ray->Direction, Depth2);
      VAddEq (IPoint, Ray->Initial);

      if (Point_In_Clip (&IPoint, Object->Clip))
	{
	push_entry(Depth2,IPoint,Object,Depth_Stack);
	Intersection_Found = TRUE;
	}
      }
    }
  return (Intersection_Found);
  }

#ifdef DB_CODE
/* Some things have been optimized. */
/*****************************************************************************
   The term VDot(Quadric->Mixed_Terms, Ray->Mixed_Initial_Initial)
   is added to the cached constant CM_Constant. Thus 3 MUL and 3 ADD
   are saved for primary rays.
   The calculation of the coefficients of the quadratic equation
   is rewritten to save 3 MUL.
 *****************************************************************************/
int Intersect_Quadric (Ray, Quadric, Depth1, Depth2)
RAY *Ray;
QUADRIC *Quadric;
DBL *Depth1, *Depth2;
{
  register DBL Square_Term, Linear_Term, Constant_Term;
  register DBL Determinant;

  Ray_Quadric_Tests++;

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

  if (Quadric->Non_Zero_Square_Term)
  {
    Square_Term = Quadric->Square_Terms.x * Ray->Direction_2.x +
		  Quadric->Square_Terms.y * Ray->Direction_2.y +
		  Quadric->Square_Terms.z * Ray->Direction_2.z +
		  Quadric->Mixed_Terms.x * Ray->Mixed_Dir_Dir.x +
		  Quadric->Mixed_Terms.y * Ray->Mixed_Dir_Dir.y +
		  Quadric->Mixed_Terms.z * Ray->Mixed_Dir_Dir.z;
  }
  else
  {
    Square_Term = 0.0;
  }

  Linear_Term = Quadric->Square_Terms.x * Ray->Initial_Direction.x +
		Quadric->Square_Terms.y * Ray->Initial_Direction.y +
		Quadric->Square_Terms.z * Ray->Initial_Direction.z +
		0.5 * (Quadric->Terms.x * Ray->Direction.x +
		       Quadric->Terms.y * Ray->Direction.y +
		       Quadric->Terms.z * Ray->Direction.z +
		       Quadric->Mixed_Terms.x * Ray->Mixed_Init_Dir.x +
		       Quadric->Mixed_Terms.y * Ray->Mixed_Init_Dir.y +
		       Quadric->Mixed_Terms.z * Ray->Mixed_Init_Dir.z);

  if (Ray == CM_Ray)
  {
    if (!Quadric->Constant_Cached)
    {
      Constant_Term = Quadric->Square_Terms.x * Ray->Initial_2.x +
		      Quadric->Square_Terms.y * Ray->Initial_2.y +
		      Quadric->Square_Terms.z * Ray->Initial_2.z +
		      Quadric->Terms.x * Ray->Initial.x +
		      Quadric->Terms.y * Ray->Initial.y +
		      Quadric->Terms.z * Ray->Initial.z +
		      Quadric->Mixed_Terms.x * Ray->Mixed_Initial_Initial.x +
		      Quadric->Mixed_Terms.y * Ray->Mixed_Initial_Initial.y +
		      Quadric->Mixed_Terms.z * Ray->Mixed_Initial_Initial.z +
		      Quadric->Constant;
      Quadric->CM_Constant = Constant_Term;
      Quadric->Constant_Cached = TRUE;
    }
    else
    {
      Constant_Term = Quadric->CM_Constant;
    }
  }
  else
  {
    Constant_Term = Quadric->Square_Terms.x * Ray->Initial_2.x +
		    Quadric->Square_Terms.y * Ray->Initial_2.y +
		    Quadric->Square_Terms.z * Ray->Initial_2.z +
		    Quadric->Terms.x * Ray->Initial.x +
		    Quadric->Terms.y * Ray->Initial.y +
		    Quadric->Terms.z * Ray->Initial.z +
		    Quadric->Mixed_Terms.x * Ray->Mixed_Initial_Initial.x +
		    Quadric->Mixed_Terms.y * Ray->Mixed_Initial_Initial.y +
		    Quadric->Mixed_Terms.z * Ray->Mixed_Initial_Initial.z +
		    Quadric->Constant;
  }

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

    Determinant = Linear_Term * Linear_Term - Square_Term * Constant_Term;

    if (Determinant < 0.0)
      return (FALSE);

    Determinant = sqrt (Determinant);

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

    *Depth1 = -2.0 * Constant_Term / 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);
}
#else
int Intersect_Quadric (Ray, Quadric, Depth1, Depth2)
RAY *Ray;
QUADRIC *Quadric;
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 (Quadric->Non_Zero_Square_Term)
    {
    VDot (Square_Term, Quadric->Square_Terms, Ray->Direction_2);
    VDot (Temp_Term, Quadric->Mixed_Terms, Ray->Mixed_Dir_Dir);
    Square_Term += Temp_Term;
    }
  else
    Square_Term = 0.0;

  VDot (Linear_Term, Quadric->Square_Terms, Ray->Initial_Direction);
  Linear_Term *= 2.0;
  VDot (Temp_Term, Quadric->Terms, Ray->Direction);
  Linear_Term += Temp_Term;
  VDot (Temp_Term, Quadric->Mixed_Terms, Ray->Mixed_Init_Dir);
  Linear_Term += Temp_Term;

  if (Ray == CM_Ray)
    if (!Quadric->Constant_Cached)
      {
      VDot (Constant_Term, Quadric->Square_Terms, Ray->Initial_2);
      VDot (Temp_Term, Quadric->Terms, Ray->Initial);
      Constant_Term +=  Temp_Term + Quadric->Constant;
      Quadric->CM_Constant = Constant_Term;
      Quadric->Constant_Cached = TRUE;
      }
    else
      Constant_Term = Quadric->CM_Constant;
  else
    {
    VDot (Constant_Term, Quadric->Square_Terms, Ray->Initial_2);
    VDot (Temp_Term, Quadric->Terms, Ray->Initial);
    Constant_Term += Temp_Term + Quadric->Constant;
    }

  VDot (Temp_Term, Quadric->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);
  }
#endif

#ifdef DB_CODE
/* Some things have been optimized. */
/*****************************************************************************
   Rewriting saves 6 MUL.
 *****************************************************************************/
int Inside_Quadric (IPoint, Object)
VECTOR *IPoint;
OBJECT *Object;
{
  QUADRIC *Quadric = (QUADRIC *) Object;
  register DBL Result;

  Result =

    IPoint->x * (Quadric->Square_Terms.x * IPoint->x +
		 Quadric->Mixed_Terms.x * IPoint->y +
		 Quadric->Terms.x) +

    IPoint->y * (Quadric->Square_Terms.y * IPoint->y +
		 Quadric->Mixed_Terms.z * IPoint->z +
		 Quadric->Terms.y) +

    IPoint->z * (Quadric->Square_Terms.z * IPoint->z +
		 Quadric->Mixed_Terms.y * IPoint->x +
		 Quadric->Terms.z) +

    Quadric->Constant;

  return (Result < Small_Tolerance);
}
#else
int Inside_Quadric (IPoint, Object)
VECTOR *IPoint;
OBJECT *Object;
  {
  QUADRIC *Quadric = (QUADRIC *) Object;
  VECTOR New_Point;
  register DBL Result, Linear_Term, Square_Term;

  VDot (Linear_Term, *IPoint, Quadric->Terms);
  Result = Linear_Term + Quadric->Constant;
  VSquareTerms (New_Point, *IPoint);
  VDot (Square_Term, New_Point, Quadric->Square_Terms);
  Result += Square_Term;
  Result += Quadric->Mixed_Terms.x * (IPoint->x) * (IPoint->y)
    + Quadric->Mixed_Terms.y * (IPoint->x) * (IPoint->z)
      + Quadric->Mixed_Terms.z * (IPoint->y) * (IPoint->z);

  if (Result < Small_Tolerance)
    return (TRUE);

  return (FALSE);
  }
#endif

#ifdef DB_CODE
/* Some things have been optimized. */
/*****************************************************************************
   Rewriting saves some memory transfers.
   Storing of 2.0 * Square_Terms would save 3 MUL.
 *****************************************************************************/
void Quadric_Normal (Result, Object, IPoint)
VECTOR *Result, *IPoint;
OBJECT *Object;
{
  QUADRIC *Quadric = (QUADRIC *) Object;
  DBL Len;

  Result->x = 2.0 * Quadric->Square_Terms.x * IPoint->x +
	      Quadric->Mixed_Terms.x * IPoint->y +
	      Quadric->Mixed_Terms.y * IPoint->z +
	      Quadric->Terms.x;

  Result->y = 2.0 * Quadric->Square_Terms.y * IPoint->y +
	      Quadric->Mixed_Terms.x * IPoint->x +
	      Quadric->Mixed_Terms.z * IPoint->z +
	      Quadric->Terms.y;

  Result->z = 2.0 * Quadric->Square_Terms.z * IPoint->z +
	      Quadric->Mixed_Terms.y * IPoint->x +
	      Quadric->Mixed_Terms.z * IPoint->y +
	      Quadric->Terms.z;

  Len = sqrt(Result->x * Result->x + Result->y * Result->y + Result->z * Result->z);

  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;
  }
}
#else
void Quadric_Normal (Result, Object, IPoint)
VECTOR *Result, *IPoint;
OBJECT *Object;
  {
  QUADRIC *Intersection_Quadric = (QUADRIC *) Object;
  VECTOR Derivative_Linear;
  DBL Len;

  VScale (Derivative_Linear, Intersection_Quadric->Square_Terms, 2.0);
  VEvaluate (*Result, Derivative_Linear, *IPoint);
  VAdd (*Result, *Result, Intersection_Quadric->Terms);

  Result->x +=
  Intersection_Quadric->Mixed_Terms.x * IPoint->y +
  Intersection_Quadric->Mixed_Terms.y * IPoint->z;


  Result->y +=
  Intersection_Quadric->Mixed_Terms.x * IPoint->x +
  Intersection_Quadric->Mixed_Terms.z * IPoint->z;

  Result->z +=
  Intersection_Quadric->Mixed_Terms.y * IPoint->x +
  Intersection_Quadric->Mixed_Terms.z * IPoint->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;
    }
  }
#endif

  void Transform_Quadric (Object, Trans)
    OBJECT *Object;
TRANSFORM *Trans;
  {
  QUADRIC *Quadric=(QUADRIC *)Object;
  MATRIX Quadric_Matrix, Transform_Transposed;

  Quadric_To_Matrix (Quadric, (MATRIX *) &Quadric_Matrix[0][0]);
  MTimes ((MATRIX *) &Quadric_Matrix[0][0], (MATRIX *) &(Trans->inverse[0][0]), (MATRIX *) &Quadric_Matrix[0][0]);
  MTranspose ((MATRIX *) &Transform_Transposed[0][0], (MATRIX *) &(Trans->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], Quadric);
#ifdef DB_CODE
  /* this has been missing. */
  recompute_bbox(&Object->Bounds, Trans);
#endif
  }

void Quadric_To_Matrix (Quadric, Matrix)
QUADRIC *Quadric;
MATRIX *Matrix;
  {
  MZero (Matrix);
  (*Matrix)[0][0] = Quadric->Square_Terms.x;
  (*Matrix)[1][1] = Quadric->Square_Terms.y;
  (*Matrix)[2][2] = Quadric->Square_Terms.z;
  (*Matrix)[0][1] = Quadric->Mixed_Terms.x;
  (*Matrix)[0][2] = Quadric->Mixed_Terms.y;
  (*Matrix)[0][3] = Quadric->Terms.x;
  (*Matrix)[1][2] = Quadric->Mixed_Terms.z;
  (*Matrix)[1][3] = Quadric->Terms.y;
  (*Matrix)[2][3] = Quadric->Terms.z;
  (*Matrix)[3][3] = Quadric->Constant;
  }

void Matrix_To_Quadric (Matrix, Quadric)
MATRIX *Matrix;
QUADRIC *Quadric;
  {
  Quadric->Square_Terms.x = (*Matrix)[0][0];
  Quadric->Square_Terms.y = (*Matrix)[1][1];
  Quadric->Square_Terms.z = (*Matrix)[2][2];
  Quadric->Mixed_Terms.x = (*Matrix)[0][1] + (*Matrix)[1][0];
  Quadric->Mixed_Terms.y = (*Matrix)[0][2] + (*Matrix)[2][0];
  Quadric->Terms.x = (*Matrix)[0][3] + (*Matrix)[3][0];
  Quadric->Mixed_Terms.z = (*Matrix)[1][2] + (*Matrix)[2][1];
  Quadric->Terms.y = (*Matrix)[1][3] + (*Matrix)[3][1];
  Quadric->Terms.z = (*Matrix)[2][3] + (*Matrix)[3][2];
  Quadric->Constant = (*Matrix)[3][3];
  }

void Translate_Quadric (Object, Vector)
OBJECT *Object;
VECTOR *Vector;
  {
  TRANSFORM Trans;

  Compute_Translation_Transform (&Trans, Vector);
  Transform_Quadric (Object, &Trans);

  }

void Rotate_Quadric (Object, Vector)
OBJECT *Object;
VECTOR *Vector;
  {
  TRANSFORM Trans;

  Compute_Rotation_Transform (&Trans, Vector);
  Transform_Quadric (Object, &Trans);
  }

void Scale_Quadric (Object, Vector)
OBJECT *Object;
VECTOR *Vector;
  {
  TRANSFORM Trans;

  Compute_Scaling_Transform (&Trans, Vector);
  Transform_Quadric (Object, &Trans);
  }

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

  VScaleEq (Quadric->Square_Terms, -1.0);
  VScaleEq (Quadric->Mixed_Terms, -1.0);
  VScaleEq (Quadric->Terms, -1.0);
  Quadric->Constant *= -1.0;
#ifdef DB_CODE
  /* Necessary for new bounding box calculation of intersections. */
  Quadric->Inverted ^= TRUE;
#endif
  }

QUADRIC *Create_Quadric()
  {
  QUADRIC *New;

  if ((New = (QUADRIC *) malloc (sizeof (QUADRIC))) == NULL)
    MAError ("quadric");

  INIT_OBJECT_FIELDS(New, QUADRIC_OBJECT, &Quadric_Methods)
    Make_Vector (&(New->Square_Terms), 1.0, 1.0, 1.0);
  Make_Vector (&(New->Mixed_Terms), 0.0, 0.0, 0.0);
  Make_Vector (&(New->Terms), 0.0, 0.0, 0.0);
  New->Constant = 1.0;
  New->CM_Constant = HUGE_VAL;
  New->Constant_Cached = FALSE;
  New->Non_Zero_Square_Term = FALSE;
#ifdef DB_CODE
  /* Necessary for new bounding box calculation of intersections. */
  New->Inverted = FALSE;
#endif

  return (New);
  }

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

  New = Create_Quadric ();
  *New = *((QUADRIC *) Object);

  return (New);
  }

void Destroy_Quadric (Object)
OBJECT *Object;
  {
  free (Object);
  }