/*****************************************************************************
*                fractal.c
*
*  This module implements the fractal sets primitive.
*
*  This file was written by Pascal Massimino.
*
*  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.
*
*****************************************************************************/

#include "frame.h"
#include "povray.h"
#include "vector.h"
#include "povproto.h"
#include "bbox.h"
#include "matrices.h"
#include "objects.h"
#include "spheres.h"
#include "fractal.h"
#include "quatern.h"
#include "hcmplx.h"



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

#ifndef Fractal_Tolerance
#define Fractal_Tolerance 1e-7
#endif



/*****************************************************************************
* Local typedefs
******************************************************************************/



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

static int All_Fractal_Intersections PARAMS((OBJECT * Object, RAY * Ray, ISTACK * Depth_Stack));
static int Inside_Fractal PARAMS((VECTOR IPoint, OBJECT * Object));
static void Fractal_Normal PARAMS((VECTOR Result, OBJECT * Object, INTERSECTION * Intersect));
static void *Copy_Fractal PARAMS((OBJECT * Object));
static void Translate_Fractal PARAMS((OBJECT * Object, VECTOR Vector, TRANSFORM *Trans));
static void Rotate_Fractal PARAMS((OBJECT * Object, VECTOR Vector, TRANSFORM *Trans));
static void Scale_Fractal PARAMS((OBJECT * Object, VECTOR Vector, TRANSFORM *Trans));
static void Transform_Fractal PARAMS((OBJECT * Object, TRANSFORM * Trans));
static void Invert_Fractal PARAMS((OBJECT * Object));
static void Destroy_Fractal PARAMS((OBJECT * Object));
static void Compute_Fractal_BBox PARAMS((FRACTAL * Fractal));

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

static METHODS Fractal_Methods =
{
  All_Fractal_Intersections,
  Inside_Fractal, Fractal_Normal,
  Copy_Fractal,
  Translate_Fractal, Rotate_Fractal,
  Scale_Fractal, Transform_Fractal, Invert_Fractal,
  Destroy_Fractal
};

static int Allocated_Iteration_Stack_Length = 0;

DBL *Sx = NULL, *Sy = NULL, *Sz = NULL, *Sw = NULL;
DBL Precision;
VECTOR Direction;

static COMPLEX_FUNCTION_METHOD Complex_Function_List[] = 
{
  /* must match STYPE list in fractal.h */
  Complex_Exp,
  Complex_Log,
  Complex_Sin,
  Complex_ASin,
  Complex_Cos,
  Complex_ACos,
  Complex_Tan,
  Complex_ATan,
  Complex_Sinh,
  Complex_ASinh,
  Complex_Cosh,
  Complex_ACosh,
  Complex_Tanh,
  Complex_ATanh,
  Complex_Pwr 
};

/*****************************************************************************
*
* FUNCTION
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   Pascal Massimino
*
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   Dec 1994 : Creation.
*
******************************************************************************/

static int All_Fractal_Intersections(Object, Ray, Depth_Stack)
OBJECT *Object;
RAY *Ray;
ISTACK *Depth_Stack;
{
  int Intersection_Found;
  int Last = 0;
  int CURRENT, NEXT;
  DBL Depth, Depth_Max;
  DBL Dist, Dist_Next, Len;

  VECTOR IPoint, Mid_Point, Next_Point, Real_Pt;
  VECTOR Real_Normal, F_Normal;
  RAY New_Ray;
  FRACTAL *Fractal = (FRACTAL *) Object;

  Increase_Counter(stats[Ray_Fractal_Tests]);

  Intersection_Found = FALSE;
  Precision = Fractal->Precision;

  /* Get into Fractal's world. */

  if (Fractal->Trans != NULL)
  {
    MInvTransDirection(Direction, Ray->Direction, Fractal->Trans);
    VDot(Len, Direction, Direction);

    if (Len == 0.0)
    {
      return (FALSE);
    }

    if (Len != 1.0)
    {
      Len = 1.0 / sqrt(Len);
      VScaleEq(Direction, Len);
    }

    Assign_Vector(New_Ray.Direction, Direction);
    MInvTransPoint(New_Ray.Initial, Ray->Initial, Fractal->Trans);
  }
  else
  {
    Assign_Vector(Direction, Ray->Direction);
    New_Ray = *Ray;
    Len = 1.0;
  }

  /* Bound fractal. */

  if (!F_Bound(&New_Ray, Fractal, &Depth, &Depth_Max))
  {
    return (FALSE);
  }

  if (Depth_Max < Fractal_Tolerance)
  {
    return (FALSE);
  }

  if (Depth < Fractal_Tolerance)
  {
    Depth = Fractal_Tolerance;
  }

  /* Jump to starting point */

  VScale(Next_Point, Direction, Depth);
  VAddEq(Next_Point, New_Ray.Initial);

  CURRENT = D_Iteration(Next_Point, Fractal, &Dist);

  /* Light ray starting inside ? */

  if (CURRENT)
  {
    VAddScaledEq(Next_Point, 2.0 * Fractal_Tolerance, Direction);

    Depth += 2.0 * Fractal_Tolerance;

    if (Depth > Depth_Max)
    {
      return (FALSE);
    }

    CURRENT = D_Iteration(Next_Point, Fractal, &Dist);
  }

  /* Ok. Trace it */

  while (Depth < Depth_Max)
  {
    /*
     * Get close to the root: Advance with Next_Point, keeping track of last
     * position in IPoint...
     */

    while (1)
    {
      if (Dist < Precision)
      {
        Dist = Precision;
      }

      Depth += Dist;

      if (Depth > Depth_Max)
      {
        if (Intersection_Found)
        {
          Increase_Counter(stats[Ray_Fractal_Tests_Succeeded]);
        }

        return (Intersection_Found);
      }

      Assign_Vector(IPoint, Next_Point);
      VAddScaledEq(Next_Point, Dist, Direction);

      NEXT = D_Iteration(Next_Point, Fractal, &Dist_Next);

      if (NEXT != CURRENT)
      {
        /* Set surface was crossed... */

        Depth -= Dist;
        break;
      }
      else
      {
        Dist = Dist_Next; /* not reached */
      }
    }

    /* then, polish the root via bisection method... */

    while (Dist > Fractal_Tolerance)
    {
      Dist *= 0.5;
      VAddScaled(Mid_Point, IPoint, Dist, Direction);

      Last = Iteration(Mid_Point, Fractal);

      if (Last == CURRENT)
      {
        Assign_Vector(IPoint, Mid_Point);

        Depth += Dist;

        if (Depth > Depth_Max)
        {
          if (Intersection_Found)
          {
            Increase_Counter(stats[Ray_Fractal_Tests_Succeeded]);
          }

          return (Intersection_Found);
        }
      }
    }

    if (CURRENT == FALSE) /* Mid_Point isn't inside the set */
    {
      VAddScaledEq(IPoint, Dist, Direction);

      Depth += Dist;

      Iteration(IPoint, Fractal);
    }
    else
    {
      if (Last != CURRENT)
      {
        Iteration(IPoint, Fractal);
      }
    }

    if (Fractal->Trans != NULL)
    {
      MTransPoint(Real_Pt, IPoint, Fractal->Trans);
      Normal_Calc(Fractal, F_Normal);
      MTransNormal(Real_Normal, F_Normal, Fractal->Trans);
    }
    else
    {
      Assign_Vector(Real_Pt, IPoint);
      Normal_Calc(Fractal, Real_Normal);
    }

    if (Point_In_Clip(Real_Pt, Object->Clip))
    {
      VNormalize(Real_Normal, Real_Normal);
      push_normal_entry(Depth * Len, Real_Pt, Real_Normal, Object, Depth_Stack);
      Intersection_Found = TRUE;

      /* If fractal isn't used with CSG we can exit now. */

      if (!(Fractal->Type & IS_CHILD_OBJECT))
      {
        break;
      }
    }

    /* Start over where work was left */

    Assign_Vector(IPoint, Next_Point);
    Dist = Dist_Next;
    CURRENT = NEXT;

  }

  if (Intersection_Found)
  {
    Increase_Counter(stats[Ray_Fractal_Tests_Succeeded]);
  }
  return (Intersection_Found);
}

/*****************************************************************************
*
* FUNCTION
*
* INPUT
*   
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   Pascal Massimino
*
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   Dec 1994 : Creation.
*
******************************************************************************/

static int Inside_Fractal(IPoint, Object)
VECTOR IPoint;
OBJECT *Object;
{
  FRACTAL *Fractal = (FRACTAL *) Object;
  int Result;
  static VECTOR New_Point;

  if (Fractal->Trans != NULL)
  {
    MInvTransPoint(New_Point, IPoint, Fractal->Trans);

    Result = Iteration(New_Point, Fractal);
  }
  else
  {
    Result = Iteration(IPoint, Fractal);
  }

  if (Fractal->Inverted)
  {
    return (!Result);
  }
  else
  {
    return (Result);
  }
}

/*****************************************************************************
*
* FUNCTION
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Pascal Massimino
*   
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   Dec 1994 : Creation.
*
******************************************************************************/

static void Fractal_Normal(Result, Object, Intersect)
OBJECT *Object;
VECTOR Result;
INTERSECTION *Intersect;
{
  Assign_Vector(Result, Intersect->INormal);
}

/*****************************************************************************
*
* FUNCTION
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Pascal Massimino
*   
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   Dec 1994 : Creation.
*
******************************************************************************/

static void Translate_Fractal(Object, Vector, Trans)
OBJECT *Object;
VECTOR Vector;
TRANSFORM *Trans;
{
  Transform_Fractal(Object, Trans);
}

/*****************************************************************************
*
* FUNCTION
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Pascal Massimino
*   
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   Dec 1994 : Creation.
*
******************************************************************************/

static void Rotate_Fractal(Object, Vector, Trans)
OBJECT *Object;
VECTOR Vector;
TRANSFORM *Trans;
{
  Transform_Fractal(Object, Trans);
}

/*****************************************************************************
*
* FUNCTION
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Pascal Massimino
*   
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   Dec 1994 : Creation.
*
******************************************************************************/

static void Scale_Fractal(Object, Vector, Trans)
OBJECT *Object;
VECTOR Vector;
TRANSFORM *Trans;
{
  Transform_Fractal(Object, Trans);
}

/*****************************************************************************
*
* FUNCTION
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   Pascal Massimino
*
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   Dec 1994 : Creation.
*   Mar 1996 : Moved call to Recompute_BBox to Compute_Fractal_BBox() (TW)
*
******************************************************************************/

static void Transform_Fractal(Object, Trans)
OBJECT *Object;
TRANSFORM *Trans;
{
  FRACTAL *Fractal = (FRACTAL *) Object;

  if (((FRACTAL *) Object)->Trans == NULL)
  {
    ((FRACTAL *) Object)->Trans = Create_Transform();
  }

  Compose_Transforms(Fractal->Trans, Trans);

  Compute_Fractal_BBox((FRACTAL *) Object);
}

/*****************************************************************************
*
* FUNCTION
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Pascal Massimino
*   
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   Dec 1994 : Creation.
*
******************************************************************************/

static void Invert_Fractal(Object)
OBJECT *Object;
{
  ((FRACTAL *) Object)->Inverted ^= TRUE;
}

/*****************************************************************************
*
* FUNCTION
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Pascal Massimino
*   
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   Dec 1994 : Creation.
*   Mar 1996 : Added call to recompute_BBox() to bottom (TW)
*
******************************************************************************/

static void Compute_Fractal_BBox(Fractal)
FRACTAL *Fractal;
{
  DBL R;

  switch (Fractal->Algebra)
  {
    case QUATERNION_TYPE:

      R = 1.0 + sqrt(Sqr(Fractal->Julia_Parm[X]) + Sqr(Fractal->Julia_Parm[Y]) + Sqr(Fractal->Julia_Parm[Z]) + Sqr(Fractal->Julia_Parm[T]));
      R += Fractal_Tolerance; /* fix bug when Julia_Parameter exactly 0 */

      if (R > 2.0)
      {
        R = 2.0;
      }

      Fractal->Exit_Value = Sqr(R);

      break;

    case HYPERCOMPLEX_TYPE:
    default:

      R = 4.0;

      Fractal->Exit_Value = 16.0;

      break;
  }

  Fractal->Radius_Squared = Sqr(R);

  Fractal->Inverted = FALSE;

  Make_BBox(Fractal->BBox, -R, -R, -R, 2.0 * R, 2.0 * R, 2.0 * R);

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

/*****************************************************************************
*
* FUNCTION
*
* INPUT
*   
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   Pascal Massimino
*
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   Dec 1994 : Creation.
*
******************************************************************************/

FRACTAL *Create_Fractal()
{
  FRACTAL *New;

  New = (FRACTAL *) POV_MALLOC(sizeof(FRACTAL), "Fractal Set");

  INIT_OBJECT_FIELDS(New, BASIC_OBJECT, &Fractal_Methods)

  New->Trans = NULL;

  Make_Vector(New->Center, 0.0, 0.0, 0.0);

  New->Julia_Parm[X] = 1.0;
  New->Julia_Parm[Y] = 0.0;
  New->Julia_Parm[Z] = 0.0;
  New->Julia_Parm[T] = 0.0;

  New->Slice[X] = 0.0;
  New->Slice[Y] = 0.0;
  New->Slice[Z] = 0.0;
  New->Slice[T] = 1.0;
  New->SliceDist = 0.0;

  New->Exit_Value = 4.0;

  New->n = 20;

  New->Precision = 1.0 / 20.0;

  New->Inverted = FALSE;

  New->Algebra = QUATERNION_TYPE;

  New->Sub_Type = SQR_STYPE;

  New->Bound = NULL;

  New->Clip = NULL;

  New->Normal_Calc_Method = NULL;
  New->Iteration_Method   = NULL;
  New->D_Iteration_Method = NULL;
  New->F_Bound_Method     = NULL;
  New->Complex_Function_Method = NULL;

  New->Radius_Squared = 0.0;
  New->exponent.x = 0.0;
  New->exponent.y = 0.0;
  
  return (New);
}

/*****************************************************************************
*
* FUNCTION
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Pascal Massimino
*   
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   Dec 1994 : Creation.
*
******************************************************************************/

static void *Copy_Fractal(Object)
OBJECT *Object;
{
  FRACTAL *New;

  New = Create_Fractal();

  *New = *((FRACTAL *) Object);

  New->Trans = Copy_Transform(((FRACTAL *) Object)->Trans);

  return (New);
}

/*****************************************************************************
*
* FUNCTION
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Pascal Massimino
*   
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   Dec 1994 : Creation.
*
******************************************************************************/

static void Destroy_Fractal(Object)
OBJECT *Object;
{
  Destroy_Transform(((FRACTAL *) Object)->Trans);
  POV_FREE(Object);
}

/*****************************************************************************
*
* FUNCTION
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Pascal Massimino
*   
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   Dec 1994 : Creation.
*
******************************************************************************/

void SetUp_Fractal(Fractal)
FRACTAL *Fractal;
{
  switch (Fractal->Algebra)
  {
    case QUATERNION_TYPE:

      switch(Fractal->Sub_Type)
      {
        case CUBE_STYPE:
          Fractal->Iteration_Method = Iteration_z3;
          Fractal->Normal_Calc_Method = Normal_Calc_z3;
          Fractal->D_Iteration_Method = D_Iteration_z3;
          break;
        case SQR_STYPE:
          Fractal->Iteration_Method = Iteration_Julia;
          Fractal->D_Iteration_Method = D_Iteration_Julia;
          Fractal->Normal_Calc_Method = Normal_Calc_Julia;
          break;
        default:
          Error("illegal function: quaternion only supports sqr and cube");
      }
      Fractal->F_Bound_Method = F_Bound_Julia;

      break;

    case HYPERCOMPLEX_TYPE:

      switch (Fractal->Sub_Type)
      {
        case RECIPROCAL_STYPE:

          Fractal->Iteration_Method = Iteration_HCompl_Reciprocal;
          Fractal->Normal_Calc_Method = Normal_Calc_HCompl_Reciprocal;
          Fractal->D_Iteration_Method = D_Iteration_HCompl_Reciprocal;
          Fractal->F_Bound_Method = F_Bound_HCompl_Reciprocal;

          break;

        case EXP_STYPE: 
        case LOG_STYPE: 
        case SIN_STYPE: 
        case ASIN_STYPE:
        case COS_STYPE: 
        case ACOS_STYPE:
        case TAN_STYPE: 
        case ATAN_STYPE:
        case SINH_STYPE:
        case ASINH_STYPE:
        case COSH_STYPE:
        case ACOSH_STYPE:
        case TANH_STYPE:
        case ATANH_STYPE:
        case PWR_STYPE: 

          Fractal->Iteration_Method = Iteration_HCompl_Func;
          Fractal->Normal_Calc_Method = Normal_Calc_HCompl_Func;
          Fractal->D_Iteration_Method = D_Iteration_HCompl_Func;
          Fractal->F_Bound_Method = F_Bound_HCompl_Func;
          Fractal->Complex_Function_Method = Complex_Function_List[Fractal->Sub_Type];

          break;

        case CUBE_STYPE:

          Fractal->Iteration_Method = Iteration_HCompl_z3;
          Fractal->Normal_Calc_Method = Normal_Calc_HCompl_z3;
          Fractal->D_Iteration_Method = D_Iteration_HCompl_z3;
          Fractal->F_Bound_Method = F_Bound_HCompl_z3;

          break;

        default:  /* SQR_STYPE or else... */

          Fractal->Iteration_Method = Iteration_HCompl;
          Fractal->Normal_Calc_Method = Normal_Calc_HCompl;
          Fractal->D_Iteration_Method = D_Iteration_HCompl;
          Fractal->F_Bound_Method = F_Bound_HCompl;

          break;
      }

      break;

    default:

      Error("Algebra unknown in fractal.");
  }

  Allocate_Iteration_Stack(Fractal->n);

  Compute_Fractal_BBox(Fractal);
}

/*****************************************************************************
*
* FUNCTION
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   Pascal Massimino
*
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   Dec 1994 : Creation.
*
******************************************************************************/

void Allocate_Iteration_Stack(n)
int n;
{
  if (n > Allocated_Iteration_Stack_Length)
  {
    Sx = (DBL *) POV_REALLOC(Sx, (n + 1) * sizeof(DBL), "x iteration stack");
    Sy = (DBL *) POV_REALLOC(Sy, (n + 1) * sizeof(DBL), "y iteration stack");
    Sz = (DBL *) POV_REALLOC(Sz, (n + 1) * sizeof(DBL), "w iteration stack");
    Sw = (DBL *) POV_REALLOC(Sw, (n + 1) * sizeof(DBL), "z iteration stack");

    Allocated_Iteration_Stack_Length = n;
  }
}

/*****************************************************************************
*
* FUNCTION
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Pascal Massimino
*   
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   Dec 1994 : Creation.
*
******************************************************************************/

void Free_Iteration_Stack()
{
  if (Sx != NULL)
  {
    POV_FREE(Sx);
    POV_FREE(Sy);
    POV_FREE(Sz);
    POV_FREE(Sw);
  }

  Sx = NULL;
  Sy = NULL;
  Sz = NULL;
  Sw = NULL;

  Allocated_Iteration_Stack_Length = 0;
}