/****************************************************************************
*                   torus.c
*
*  This module implements functions that manipulate torii.
*
*  This module was written by Dieter Bayer [DB].
*
*  from Persistence of Vision(tm) Ray Tracer
*  Copyright 1996 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.
*
*****************************************************************************/

/****************************************************************************
*
*  Explanation:
*
*  ---
*
*  June 1994 : Creation. [DB]
*
*****************************************************************************/

#include "frame.h"
#include "povray.h"
#include "vector.h"
#include "povproto.h"
#include "bbox.h"
#include "polysolv.h"
#include "matrices.h"
#include "objects.h"
#include "torus.h"



/*****************************************************************************
* Local preprocessor defines
******************************************************************************/

/* Minimal depth for a valid intersection. */

#define DEPTH_TOLERANCE 1.0e-4

/* Tolerance used for order reduction during root finding. */

#define ROOT_TOLERANCE 1.0e-4



/*****************************************************************************
* Static functions
******************************************************************************/

static int intersect_torus PARAMS((RAY *Ray, TORUS *Torus, DBL *Depth));
static int   All_Torus_Intersections PARAMS((OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack));
static int   Inside_Torus PARAMS((VECTOR point, OBJECT *Object));
static void  Torus_Normal PARAMS((VECTOR Result, OBJECT *Object, INTERSECTION *Inter));
static void  *Copy_Torus PARAMS((OBJECT *Object));
static void  Translate_Torus PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
static void  Rotate_Torus PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
static void  Scale_Torus PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
static void  Transform_Torus PARAMS((OBJECT *Object, TRANSFORM *Trans));
static void  Invert_Torus PARAMS((OBJECT *Object));
static void  Destroy_Torus PARAMS((OBJECT *Object));



/*****************************************************************************
* Local variables
******************************************************************************/

static METHODS Torus_Methods =
{
  All_Torus_Intersections, Inside_Torus, Torus_Normal,
  Copy_Torus, Translate_Torus, Rotate_Torus,
  Scale_Torus, Transform_Torus, Invert_Torus, Destroy_Torus
};


/*****************************************************************************
*
* FUNCTION
*
*   All_Torus_Intersections
*
* INPUT
*
*   Object      - Object
*   Ray         - Ray
*   Depth_Stack - Intersection stack
*   
* OUTPUT
*
*   Depth_Stack
*   
* RETURNS
*
*   int - TRUE, if an intersection was found
*   
* AUTHOR
*
*   Dieter Bayer
*   
* DESCRIPTION
*
*   Determine ray/torus intersection and clip intersection found.
*
* CHANGES
*
*   Jun 1994 : Creation.
*
******************************************************************************/

static int All_Torus_Intersections(Object, Ray, Depth_Stack)
OBJECT *Object;
RAY *Ray;
ISTACK *Depth_Stack;
{
  int i, max_i, Found;
  DBL Depth[4];
  VECTOR IPoint;

  Found = FALSE;

  if ((max_i = intersect_torus(Ray, (TORUS *)Object, Depth)) > 0)
  {
    for (i = 0; i < max_i; i++)
    {
      if ((Depth[i] > DEPTH_TOLERANCE) && (Depth[i] < Max_Distance))
      {
        VEvaluateRay(IPoint, Ray->Initial, Depth[i], Ray->Direction);

        if (Point_In_Clip(IPoint, Object->Clip))
        {
          push_entry(Depth[i], IPoint, Object, Depth_Stack);

          Found = TRUE;
        }
      }
    }
  }

  return(Found);
}



/*****************************************************************************
*
* FUNCTION
*
*   intersect_torus
*
* INPUT
*
*   Ray   - Ray
*   Torus - Torus
*   Depth - Intersections found
*   
* OUTPUT
*
*   Depth
*   
* RETURNS
*
*   int - Number of intersections found
*   
* AUTHOR
*
*   Dieter Bayer
*   
* DESCRIPTION
*
*   Determine ray/torus intersection.
*
*   Note that the torus is rotated about the y-axis!
*
* CHANGES
*
*   Jun 1994 : Creation.
*
******************************************************************************/

static int intersect_torus(Ray, Torus, Depth)
RAY *Ray;
TORUS *Torus;
DBL *Depth;
{
  int i, n;
  DBL len, R2, Py2, Dy2, PDy2, k1, k2;
  DBL y1, y2, r1, r2;
  DBL c[5], r[4];
  VECTOR P, D;

  Increase_Counter(stats[Ray_Torus_Tests]);

  /* Transform the ray into the torus space. */

  MInvTransPoint(P, Ray->Initial, Torus->Trans);

  MInvTransDirection(D, Ray->Direction, Torus->Trans);

  VLength(len, D);

  VInverseScaleEq(D, len);

  i = 0;

  y1 = -Torus->r;
  y2 =  Torus->r;
  r1 = Sqr(Torus->R - Torus->r);
  r2 = Sqr(Torus->R + Torus->r);

#ifdef TORUS_EXTRA_STATS
  Increase_Counter(stats[Torus_Bound_Tests]);
#endif

  if (Test_Thick_Cylinder(P, D, y1, y2, r1, r2))
  {
#ifdef TORUS_EXTRA_STATS
    Increase_Counter(stats[Torus_Bound_Tests_Succeeded]);
#endif

    R2   = Sqr(Torus->R);
    r2   = Sqr(Torus->r);

    Py2  = P[Y] * P[Y];
    Dy2  = D[Y] * D[Y];
    PDy2 = P[Y] * D[Y];

    k1   = P[X] * P[X] + P[Z] * P[Z] + Py2 - R2 - r2;
    k2   = P[X] * D[X] + P[Z] * D[Z] + PDy2;

    c[0] = 1.0;

    c[1] = 4.0 * k2;

    c[2] = 2.0 * (k1 + 2.0 * (k2 * k2 + R2 * Dy2));

    c[3] = 4.0 * (k2 * k1 + 2.0 * R2 * PDy2);

    c[4] = k1 * k1 + 4.0 * R2 * (Py2 - r2);

    n = Solve_Polynomial(4, c, r, Test_Flag(Torus, STURM_FLAG), ROOT_TOLERANCE);

    while(n--)
    {
      Depth[i++] = r[n] / len;
    }
  }

  if (i)
  {
    Increase_Counter(stats[Ray_Torus_Tests_Succeeded]);
  }

  return(i);
}



/*****************************************************************************
*
* FUNCTION
*
*   Inside_Torus
*
* INPUT
*
*   IPoint - Intersection point
*   Object - Object
*   
* OUTPUT
*   
* RETURNS
*
*   int - TRUE if inside
*   
* AUTHOR
*
*   Dieter Bayer
*   
* DESCRIPTION
*
*   Test if a point lies inside the torus.
*
* CHANGES
*
*   Jun 1994 : Creation.
*
******************************************************************************/

static int Inside_Torus(IPoint, Object)
VECTOR IPoint;
OBJECT *Object;
{
  DBL r, r2;
  VECTOR P;
  TORUS *Torus = (TORUS *)Object;

  /* Transform the point into the torus space. */

  MInvTransPoint(P, IPoint, Torus->Trans);

  r  = sqrt(Sqr(P[X]) + Sqr(P[Z]));

  r2 = Sqr(P[Y]) + Sqr(r - Torus->R);

  if (r2 <= Sqr(Torus->r))
  {
    return(!Test_Flag(Torus, INVERTED_FLAG));
  }
  else
  {
    return(Test_Flag(Torus, INVERTED_FLAG));
  }
}



/*****************************************************************************
*
* FUNCTION
*
*   Torus_Normal
*
* INPUT
*
*   Result - Normal vector
*   Object - Object
*   Inter  - Intersection found
*   
* OUTPUT
*
*   Result
*   
* RETURNS
*   
* AUTHOR
*
*   Dieter Bayer
*   
* DESCRIPTION
*
*   Calculate the normal of the torus in a given point.
*
* CHANGES
*
*   Jun 1994 : Creation.
*
******************************************************************************/

static void Torus_Normal(Result, Object, Inter)
OBJECT *Object;
VECTOR Result;
INTERSECTION *Inter;
{
  DBL dist;
  VECTOR P, N, M;
  TORUS *Torus = (TORUS *)Object;

  /* Transform the point into the torus space. */

  MInvTransPoint(P, Inter->IPoint, Torus->Trans);

  /* Get normal from derivatives. */

  dist = sqrt(P[X] * P[X] + P[Z] * P[Z]);

  if (dist > EPSILON)
  {
    M[X] = Torus->R * P[X] / dist;
    M[Y] = 0.0;
    M[Z] = Torus->R * P[Z] / dist;
  }
  else
  {
    Make_Vector(M, 0.0, 0.0, 0.0);
  }

  VSub(N, P, M);

  /* Transform the normalt out of the torus space. */

  MTransNormal(Result, N, Torus->Trans);

  VNormalize(Result, Result);
}



/*****************************************************************************
*
* FUNCTION
*
*   Translate_Torus
*
* INPUT
*
*   Object - Object
*   Vector - Translation vector
*   
* OUTPUT
*
*   Object
*   
* RETURNS
*   
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   Translate a torus.
*
* CHANGES
*
*   Jun 1994 : Creation.
*
******************************************************************************/

static void Translate_Torus(Object, Vector, Trans)
OBJECT *Object;
VECTOR Vector;
TRANSFORM *Trans;
{
  Transform_Torus(Object, Trans);
}



/*****************************************************************************
*
* FUNCTION
*
*   Rotate_Torus
*
* INPUT
*
*   Object - Object
*   Vector - Rotation vector
*   
* OUTPUT
*
*   Object
*   
* RETURNS
*   
* AUTHOR
*
*   Dieter Bayer
*   
* DESCRIPTION
*
*   Rotate a torus.
*
* CHANGES
*
*   Jun 1994 : Creation.
*
******************************************************************************/

static void Rotate_Torus(Object, Vector, Trans)
OBJECT *Object;
VECTOR Vector;
TRANSFORM *Trans;
{
  Transform_Torus(Object, Trans);
}



/*****************************************************************************
*
* FUNCTION
*
*   Scale_Torus
*
* INPUT
*
*   Object - Object
*   Vector - Scaling vector
*   
* OUTPUT
*
*   Object
*   
* RETURNS
*   
* AUTHOR
*
*   Dieter Bayer
*   
* DESCRIPTION
*
*   Scale a torus.
*
* CHANGES
*
*   Jun 1994 : Creation.
*
******************************************************************************/

static void Scale_Torus(Object, Vector, Trans)
OBJECT *Object;
VECTOR Vector;
TRANSFORM *Trans;
{
  Transform_Torus(Object, Trans);
}



/*****************************************************************************
*
* FUNCTION
*
*   Transform_Torus
*
* INPUT
*
*   Object - Object
*   Trans  - Transformation to apply
*   
* OUTPUT
*
*   Object
*   
* RETURNS
*   
* AUTHOR
*
*   Dieter Bayer
*   
* DESCRIPTION
*
*   Transform a torus and recalculate its bounding box.
*
* CHANGES
*
*   Jun 1994 : Creation.
*
******************************************************************************/

static void Transform_Torus(Object, Trans)
OBJECT *Object;
TRANSFORM *Trans;
{
  Compose_Transforms(((TORUS *)Object)->Trans, Trans);

  Compute_Torus_BBox((TORUS *)Object);
}



/*****************************************************************************
*
* FUNCTION
*
*   Invert_Torus
*
* INPUT
*
*   Object - Object
*   
* OUTPUT
*
*   Object
*   
* RETURNS
*   
* AUTHOR
*
*   Dieter Bayer
*   
* DESCRIPTION
*
*   Invert a torus.
*
* CHANGES
*
*   Jun 1994 : Creation.
*
******************************************************************************/

static void Invert_Torus(Object)
OBJECT *Object;
{
  Invert_Flag(Object, INVERTED_FLAG);
}



/*****************************************************************************
*
* FUNCTION
*
*   Create_Torus
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*
*   TORUS * - new torus
*   
* AUTHOR
*
*   Dieter Bayer
*   
* DESCRIPTION
*
*   Create a new torus.
*
* CHANGES
*
*   Jun 1994 : Creation.
*
******************************************************************************/

TORUS *Create_Torus()
{
  TORUS *New;

  New = (TORUS *)POV_MALLOC(sizeof(TORUS), "torus");

  INIT_OBJECT_FIELDS(New,TORUS_OBJECT,&Torus_Methods)

  New->Trans = Create_Transform();

  New->R =
  New->r = 0.0;

  return (New);
}



/*****************************************************************************
*
* FUNCTION
*
*   Copy_Torus
*
* INPUT
*
*   Object - Object
*   
* OUTPUT
*   
* RETURNS
*
*   void * - New torus
*   
* AUTHOR
*
*   Dieter Bayer
*   
* DESCRIPTION
*
*   Copy a torus.
*
* CHANGES
*
*   Jun 1994 : Creation.
*
*   Sep 1994 : fixed memory leakage [DB]
*
******************************************************************************/

static void *Copy_Torus(Object)
OBJECT *Object;
{
  TORUS *New, *Torus = (TORUS *)Object;

  New = Create_Torus();

  /* Get rid of the transformation created in Create_Torus(). */

  Destroy_Transform(New->Trans);

  /* Copy torus. */

  *New = *Torus;

  New->Trans = Copy_Transform(Torus->Trans);

  return (New);
}



/*****************************************************************************
*
* FUNCTION
*
*   Destroy_Torus
*
* INPUT
*
*   Object - Object
*   
* OUTPUT
*
*   Object
*   
* RETURNS
*   
* AUTHOR
*
*   Dieter Bayer
*   
* DESCRIPTION
*
*   Destroy a torus.
*
* CHANGES
*
*   Jun 1994 : Creation.
*
******************************************************************************/

static void Destroy_Torus (Object)
OBJECT *Object;
{
  Destroy_Transform(((TORUS *)Object)->Trans);

  POV_FREE (Object);
}



/*****************************************************************************
*
* FUNCTION
*
*   Compute_Torus_BBox
*
* INPUT
*
*   Torus - Torus
*   
* OUTPUT
*
*   Torus
*   
* RETURNS
*   
* AUTHOR
*
*   Dieter Bayer
*   
* DESCRIPTION
*
*   Calculate the bounding box of a torus.
*
* CHANGES
*
*   Jun 1994 : Creation.
*
******************************************************************************/

void Compute_Torus_BBox(Torus)
TORUS *Torus;
{
  DBL r1, r2;

  r1 = Torus->r;
  r2 = Torus->R + Torus->r;

  Make_BBox(Torus->BBox, -r2, -r1, -r2, 2.0 * r2, 2.0 * r1, 2.0 * r2);

  Recompute_BBox(&Torus->BBox, Torus->Trans);
}



/*****************************************************************************
*
* FUNCTION
*
*   Test_Thick_Cylinder
*
* INPUT
*
*   P  - Ray initial
*   D  - Ray direction
*   h1 - Height 1
*   h2 - Height 2
*   r1 - Square of inner radius
*   r2 - Square of outer radius
*   
* OUTPUT
*   
* RETURNS
*
*   int - TRUE, if hit
*   
* AUTHOR
*
*   Dieter Bayer
*   
* DESCRIPTION
*
*   Test if a given ray defined in the lathe's coordinate system
*   intersects a "thick" cylinder (rotated about y-axis).
*
* CHANGES
*
*   Jun 1994 : Creation.
*
******************************************************************************/

int Test_Thick_Cylinder(P, D, h1, h2, r1, r2)
VECTOR P, D;
DBL h1, h2, r1, r2;
{
  DBL a, b, c, d;
  DBL u, v, k, r, h;

  if (fabs(D[Y]) < EPSILON)
  {
    if ((P[Y] < h1) || (P[Y] > h2))
    {
      return(FALSE);
    }
  }
  else
  {
    /* Intersect ray with the cap-plane. */

    k = (h2 - P[Y]) / D[Y];

    u = P[X] + k * D[X];
    v = P[Z] + k * D[Z];

    if ((k > EPSILON) && (k < Max_Distance))
    {
      r = u * u + v * v;

      if ((r >= r1) && (r <= r2))
      {
        return(TRUE);
      }
    }

    /* Intersectionersect ray with the base-plane. */

    k = (h1 - P[Y]) / D[Y];

    u = P[X] + k * D[X];
    v = P[Z] + k * D[Z];

    if ((k > EPSILON) && (k < Max_Distance))
    {
      r = u * u + v * v;

      if ((r >= r1) && (r <= r2))
      {
        return(TRUE);
      }
    }
  }

  a = D[X] * D[X] + D[Z] * D[Z];

  if (a > EPSILON)
  {
    /* Intersect with outer cylinder. */

    b = P[X] * D[X] + P[Z] * D[Z];

    c = P[X] * P[X] + P[Z] * P[Z] - r2;

    d = b * b - a * c;

    if (d >= 0.0)
    {
      d = sqrt(d);

      k = (-b + d) / a;

      if ((k > EPSILON) && (k < Max_Distance))
      {
        h = P[Y] + k * D[Y];

        if ((h >= h1) && (h <= h2))
        {
          return(TRUE);
        }
      }

      k = (-b - d) / a;

      if ((k > EPSILON) && (k < Max_Distance))
      {
        h = P[Y] + k * D[Y];

        if ((h >= h1) && (h <= h2))
        {
          return(TRUE);
        }
      }
    }

    /* Intersect with inner cylinder. */

    c = P[X] * P[X] + P[Z] * P[Z] - r1;

    d = b * b - a * c;

    if (d >= 0.0)
    {
      d = sqrt(d);

      k = (-b + d) / a;

      if ((k > EPSILON) && (k < Max_Distance))
      {
        h = P[Y] + k * D[Y];

        if ((h >= h1) && (h <= h2))
        {
          return(TRUE);
        }
      }

      k = (-b - d) / a;

      if ((k > EPSILON) && (k < Max_Distance))
      {
        h = P[Y] + k * D[Y];

        if ((h >= h1) && (h <= h2))
        {
          return(TRUE);
        }
      }
    }
  }

  return(FALSE);
}