/****************************************************************************
*                blob.c
*		  
*  This module contains the code for the blob shape.
*
*  This file was written by Alexander Enzmann.	He wrote the code for
*  blobs and generously provided us these enhancements.
*
*  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"

#ifndef min
#define min(x,y) ((x)<(y)?(x):(y))
#endif
#ifndef max
#define max(x,y) ((x)>(y)?(x):(y))
#endif

METHODS Blob_Methods =
  { 
  All_Blob_Intersections,
  Inside_Blob, Blob_Normal,
  Copy_Blob,
  Translate_Blob, Rotate_Blob, Scale_Blob, Transform_Blob,
  Invert_Blob, Destroy_Blob
};

extern long Ray_Blob_Tests, Ray_Blob_Tests_Succeeded;
extern int Shadow_Test_Flag;

static int determine_influences PARAMS((VECTOR *P, VECTOR *D, BLOB *blob, DBL mindist));
static DBL calculate_field_value PARAMS ((OBJECT *obj, VECTOR *Pos));

#ifndef Blob_Tolerance 
#define Blob_Tolerance 1.0e-3
#endif

#define COEFF_LIMIT 1.0e-20
#define INSIDE_TOLERANCE 1.0e-6

/* Starting with the density function: (1-r^2)^2, we have a field
   that varies in strength from 1 at r = 0 to 0 at r = 1.  By
   substituting r/rad for r, we can adjust the range of influence
   of a particular component.  By multiplication by coeff, we can
   adjust the amount of total contribution, giving the formula:
      coeff * (1 - (r/rad)^2)^2
   This varies in strength from coeff at r = 0, to 0 at r = rad. */
void
MakeBlob(blob, threshold, bloblist, npoints, sflag)
BLOB *blob;
DBL threshold;
blobstackptr bloblist;
int npoints;
int sflag;
  {
  unsigned i;
  DBL rad, coeff;
  blobstackptr temp;
  VECTOR mins, maxs;

  if (npoints < 1) 
    Error("Need at least one component in a blob.");

  blob->threshold = threshold;
  blob->list = (Blob_Element **)malloc(npoints*sizeof(Blob_Element *));
  if (blob->list == NULL) 
    MAError("blob data");
  for (i=0;i < (unsigned)npoints;i++) 
    {
    blob->list[i] = (Blob_Element *)malloc(sizeof(Blob_Element));
    if (blob->list[i] == NULL)
      MAError("blob data");
    }

  blob->count = npoints;
  blob->Sturm_Flag = sflag;

  /* Initialize the blob data */
  for(i=0;i < (unsigned)npoints;i++) 
    {
    temp = bloblist;
    if (fabs(temp->elem.coeffs[2]) < EPSILON ||
      temp->elem.radius2 < EPSILON) 
      {
      perror("Degenerate blob element\n");
      }
    /* Store blob specific information */
    rad = temp->elem.radius2;
    rad *= rad;
    coeff = temp->elem.coeffs[2];
    blob->list[i]->radius2   = rad;
    blob->list[i]->coeffs[2] = coeff;
    blob->list[i]->coeffs[1] = -(2.0 * coeff) / rad;
    blob->list[i]->coeffs[0] = coeff / (rad * rad);
    blob->list[i]->pos.x = temp->elem.pos.x;
    blob->list[i]->pos.y = temp->elem.pos.y;
    blob->list[i]->pos.z = temp->elem.pos.z;

    rad = temp->elem.radius2;
    if (i == 0) 
      {
      /* First component, just set the bounds */
      Make_Vector(&mins,
        temp->elem.pos.x-rad,
        temp->elem.pos.y-rad,
        temp->elem.pos.z-rad)
        Make_Vector(&maxs,
          temp->elem.pos.x+rad,
          temp->elem.pos.y+rad,
          temp->elem.pos.z+rad)
          }
    else 
      {
      /* Check min/max on the bounds */
        mins.x = min(mins.x, temp->elem.pos.x - rad);
      mins.y = min(mins.y, temp->elem.pos.y - rad);
      mins.z = min(mins.z, temp->elem.pos.z - rad);
      maxs.x = max(maxs.x, temp->elem.pos.x + rad);
      maxs.y = max(maxs.y, temp->elem.pos.y + rad);
      maxs.z = max(maxs.z, temp->elem.pos.z + rad);
      }

    bloblist = bloblist->next;
    free(temp);
    }

  blob->Bounds.Lower_Left = mins;
  VSub(blob->Bounds.Lengths, maxs, mins);

  /*  Allocate memory for intersection intervals */
  npoints *= 2;
  blob->intervals = (Blob_Interval *)malloc(npoints*sizeof(Blob_Interval));
  if (blob->intervals == NULL) 
    MAError("blob data");
  }

/* Make a sorted list of points along the ray that the various blob
   components start and stop adding their influence.  It would take
   a very complex blob (with many components along the current ray)
   to warrant the overhead of using a faster sort technique. */
static int
determine_influences(P, D, blob, mindist)
VECTOR *P, *D;
BLOB *blob;
DBL mindist;
  {
  int i, j, k, cnt;
  DBL b, t, t0, t1, disc;
  VECTOR V;
  Blob_Interval *intervals = blob->intervals;

  cnt = 0;
  for (i=0;i<blob->count;i++) 
    {
    /* Use standard sphere intersection routine
         to determine where the ray hits the volume
         of influence of each component of the blob. */
    VSub(V, blob->list[i]->pos, *P);
    VDot(b, V, *D);
    VDot(t, V, V);
    disc = b * b - t + blob->list[i]->radius2;
    if (disc < EPSILON)
      continue;
    disc = sqrt(disc);
    t1 = b + disc;
    if (t1 < mindist) t1 = 0.0;
    t0 = b - disc;
    if (t0 < mindist) t0 = 0.0;
    if (t1 == t0) continue;
    else if (t1 < t0) 
      {
      disc = t0;
      t0 = t1;
      t1 = disc;
      }

    /* Store the points of intersection of this
         blob with the ray.  Keep track of: whether
         this is the start or end point of the hit,
         which component was pierced by the ray,
         and the point along the ray that the
         hit occured at. */
    for (k=0;k<cnt && t0 > intervals[k].bound;k++);
    if (k<cnt) 
      {
      /* This hit point is smaller than one that
            already exists - bump the rest and insert
            it here */
      for (j=cnt;j>k;j--)
        memcpy(&intervals[j], &intervals[j-1],
          sizeof(Blob_Interval));
      intervals[k].type  = 0;
      intervals[k].index = i;
      intervals[k].bound = t0;
      cnt++;
      for (k=k+1;k<cnt && t1 > intervals[k].bound;k++);
      if (k<cnt) 
        {
        for (j=cnt;j>k;j--)
          memcpy(&intervals[j], &intervals[j-1],
            sizeof(Blob_Interval));
        intervals[k].type  = 1;
        intervals[k].index = i;
        intervals[k].bound = t1;
        }
      else 
        {
        intervals[cnt].type  = 1;
        intervals[cnt].index = i;
        intervals[cnt].bound = t1;
        }
      cnt++;
      }
    else 
      {
      /* Just plop the start and end points at
            the end of the list */
      intervals[cnt].type  = 0;
      intervals[cnt].index = i;
      intervals[cnt].bound = t0;
      cnt++;
      intervals[cnt].type  = 1;
      intervals[cnt].index = i;
      intervals[cnt].bound = t1;
      cnt++;
      }
    }
  return cnt;
  }

  /* Calculate the field value of a blob - the position vector
   "Pos" must already have been transformed into blob space. */
  static DBL
  calculate_field_value(obj, Pos)
    OBJECT *obj;
VECTOR *Pos;
  {
  int i;
  DBL len, density;
  VECTOR V;
  Blob_Element *ptr;
  BLOB *blob = (BLOB *)obj;

  density = 0.0;
  for (i=0;i<blob->count;i++) 
    {
    ptr = blob->list[i];
    VSub(V, ptr->pos, *Pos);
    VDot(len, V, V);
    if (len < ptr->radius2) 
      {
      /* Inside the radius of influence of this
            component, add it's contribution */
      density += len * (len * ptr->coeffs[0] +
        ptr->coeffs[1]) +
      ptr->coeffs[2];
      }
    }
  return density;
  }

/* Generate intervals of influence of each component.  After these
   are made, determine their aggregate effect on the ray.  As the
   individual intervals are checked, a quartic is generated
   that represents the density at a particular point on the ray.

   After making the substitutions in MakeBlob, there is a formula
   for each component that has the form:
   
      c0 * r^4 + c1 * r^2 + c2.
   
   In order to determine the influence on the ray of all of the
   individual components, we start by determining the distance
   from any point on the ray to the specified point.  This can
   be found using the pythagorean theorem, using C as the center
   of this component, P as the start of the ray, and D as the
   direction of travel of the ray:

      r^2 = (t * D + P - C) . (t * D + P - C)

   we insert this equation for each appearance of r^2 in the
   components' formula, giving:

      r^2 = D.D t^2 + 2 t D . (P - C) + (P - C) . (P - C)

   Since the direction vector has been normalized, D.D = 1.
   Using the substitutions:

      t0 = (P - C) . (P - C),
      t1 = D . (P - C)

   We can write the formula as:

      r^2 = t0 + 2 t t1 + t^2

   Taking r^2 and substituting into the formula for this component
   of the blob we get the formula:

      density = c0 * (r^2)^2 + c1 * r^2 + c2,

   or:

      density = c0 * (t0 + 2 t t1 + t^2)^2 +
                c1 * (t0 + 2 t t1 + t^2) +
                c2

   Expanding terms and collecting with respect to "t" gives:
      t^4 * c0 +
      t^3 * 4 c0 t1 +
      t^2 * (c1 + 2 * c0 t0 + 4 c0 t1^2)
      t   * 2 (c1 t1 + 2 c0 t0 t1) +
            c2 + c1*t0 + c0*t0^2

   This formula can now be solved for "t" by any of the quartic
   root solvers that are available.
*/
int All_Blob_Intersections(Object, Ray, Depth_Stack)
OBJECT *Object;
RAY *Ray;
ISTACK *Depth_Stack;
  {
  BLOB *blob = (BLOB *)Object;
  DBL dist, len, *tcoeffs, coeffs[5], roots[4];
  int i, j, cnt;
  VECTOR P, D, V;
  int root_count, in_flag;
  Blob_Element *element;
  DBL t0, t1, c0, c1, c2;
  VECTOR IPoint, dv;
  Blob_Interval *intervals = blob->intervals;
  int Intersection_Found = FALSE;

  Ray_Blob_Tests++;

  /* Transform the ray into the blob space */
  if (blob->Trans != NULL) 
    {
    MInvTransPoint(&P, &Ray->Initial, blob->Trans);
    MInvTransDirection(&D, &Ray->Direction, blob->Trans);
    }
  else 
    {
    P = Ray->Initial;
    D = Ray->Direction;
    }

    len = sqrt(D.x * D.x + D.y * D.y + D.z * D.z);
  if (len == 0.0)
    return 0;
  else 
    {
    D.x /= len;
    D.y /= len;
    D.z /= len;
    }

    /* Figure out the intervals along the ray where each
   component of the blob has an effect. */
    if ((cnt = determine_influences(&P, &D, blob, 0.01)) == 0)
      /* Ray doesn't hit the sphere of influence of any of
      its component elements */
      return 0;

  /* Clear out the coefficients */
  for (i=0;i<4;i++) coeffs[i] = 0.0;
  coeffs[4] = -blob->threshold;

  /* Step through the list of influence points, adding the
      influence of each blob component as it appears */
  for (i=0,in_flag=0;i<cnt;i++) 
    {
    if (intervals[i].type == 0) 
      {
      /* Something is just starting to influence the ray,
            so calculate its coefficients and add them
         into the pot. */
      in_flag++;
      element = blob->list[intervals[i].index];

      VSub(V, P, element->pos);
      c0 = element->coeffs[0];
      c1 = element->coeffs[1];
      c2 = element->coeffs[2];
      VDot(t0, V, V);
      VDot(t1, V, D);
      tcoeffs = &(element->tcoeffs[0]);

      tcoeffs[0] = c0;
      tcoeffs[1] = 4.0 * c0 * t1;
      tcoeffs[2] = 2.0 * c0 * (2.0 * t1 * t1 + t0) + c1;
      tcoeffs[3] = 2.0 * t1 * (2.0 * c0 * t0 + c1);
      tcoeffs[4] = c0 * t0 * t0 + c1 * t0 + c2;

      for (j=0;j<5;j++) coeffs[j] += tcoeffs[j];
      }
    else 
      {
      /* We are losing the influence of a component, so
            subtract off its coefficients */
      tcoeffs = &(blob->list[intervals[i].index]->tcoeffs[0]);
      for (j=0;j<5;j++) coeffs[j] -= tcoeffs[j];
      if (--in_flag == 0)
        /* None of the components are currently affecting
               the ray - skip ahead. */
        continue;
      }

    /* Figure out which root solver to use */
    if (blob->Sturm_Flag == 0)
      /* Use Ferrari's method */
      root_count = solve_quartic(coeffs, &roots[0]);
    else
      /* Sturm sequences */
      if (fabs(coeffs[0]) < COEFF_LIMIT)
        if (fabs(coeffs[1]) < COEFF_LIMIT)
          root_count = solve_quadratic(&coeffs[2], &roots[0]);
        else
          root_count = polysolve(3, &coeffs[1], &roots[0]);
      else
        root_count = polysolve(4, coeffs, &roots[0]);

    /* See if any of the roots are valid */
    for(j=0;j<root_count;j++) 
      {
      dist = roots[j];
      /* First see if the root is in the interval of influence of
            the currently active components of the blob */
      if ((dist >= intervals[i].bound) &&
        (dist <= intervals[i+1].bound) &&
        (dist > Blob_Tolerance)) 
        {
        VScale(IPoint, D, dist);
        VAdd(IPoint, IPoint, P);
        /* Transform the point into world space */
        if (blob->Trans != NULL)
          MTransPoint(&IPoint, &IPoint, blob->Trans);
        VSub(dv, IPoint, Ray->Initial);
        VLength(len, dv);
        if (Point_In_Clip(&IPoint, Object->Clip)) 
          {
          push_entry(len,IPoint,Object,Depth_Stack);
          Intersection_Found = TRUE;
          }
        }
      }
    }
  if(Intersection_Found)
    Ray_Blob_Tests_Succeeded++;
  return Intersection_Found;
  }

/* Calculate the density at this point, then compare to
   the threshold to see if we are in or out of the blob */
int Inside_Blob (Test_Point, Object)
VECTOR *Test_Point;
OBJECT *Object;
  {
  VECTOR New_Point;
  BLOB *blob = (BLOB *) Object;

  /* Transform the point into blob space */
  if (blob->Trans != NULL)
    MInvTransPoint(&New_Point, Test_Point, blob->Trans);
  else
    New_Point = *Test_Point;

  if (calculate_field_value(Object, &New_Point) >
    blob->threshold-INSIDE_TOLERANCE)
    return ((int) 1-blob->Inverted);
  else
    return ((int) blob->Inverted);
  }

void Blob_Normal (Result, Object, IPoint)
OBJECT *Object;
VECTOR *Result, *IPoint;
  {
  VECTOR New_Point, V;
  int i;
  DBL dist, val;
  BLOB *blob = (BLOB *) Object;
  Blob_Element *temp;

  /* Transform the point into the blobs space */
  if (blob->Trans != NULL)
    MInvTransPoint(&New_Point, IPoint, blob->Trans);
  else
    New_Point = *IPoint;

  Make_Vector(Result, 0, 0, 0);

  /* For each component that contributes to this point, add
      its bit to the normal */
  for(i=0;i<blob->count;i++) 
    {
    temp = blob->list[i];
    V.x = New_Point.x - temp->pos.x;
    V.y = New_Point.y - temp->pos.y;
    V.z = New_Point.z - temp->pos.z;
    dist = (V.x * V.x + V.y * V.y + V.z * V.z);

    if (dist <= temp->radius2) 
      {
      val = -2.0 * (2.0 * temp->coeffs[0] * dist +
        temp->coeffs[1]);
      Result->x += val * V.x;
      Result->y += val * V.y;
      Result->z += val * V.z;
      }
    }
  val = (Result->x * Result->x + Result->y * Result->y +
    Result->z * Result->z);
  if (val < EPSILON) 
    {
    Result->x = 1.0;
    Result->y = 0.0;
    Result->z = 0.0;
    }
  else 
    {
    val = 1.0 / sqrt(val);
    Result->x *= val;
    Result->y *= val;
    Result->z *= val;
    }

    /* Transform back to world space */
    if (blob->Trans != NULL)
      MTransNormal(Result, Result, blob->Trans);
  VNormalize(*Result, *Result);
  }

void Translate_Blob (Object, Vector)
OBJECT *Object;
VECTOR *Vector;
  {
  TRANSFORM Trans;
  Compute_Translation_Transform(&Trans, Vector);
  Transform_Blob(Object, &Trans);
  }

void Rotate_Blob (Object, Vector)
OBJECT *Object;
VECTOR *Vector;
  {
  TRANSFORM Trans;
  Compute_Rotation_Transform(&Trans, Vector);
  Transform_Blob(Object, &Trans);
  }

void Scale_Blob (Object, Vector)
OBJECT *Object;
VECTOR *Vector;
  {
  TRANSFORM Trans;
  Compute_Scaling_Transform(&Trans, Vector);
  Transform_Blob(Object, &Trans);
  }

void Transform_Blob(Object,Trans)
OBJECT *Object;
TRANSFORM *Trans;
  {
  if (((BLOB *)Object)->Trans == NULL)
    ((BLOB *)Object)->Trans = Create_Transform();
  recompute_bbox(&Object->Bounds, Trans);
  Compose_Transforms(((BLOB *)Object)->Trans, Trans);
  }

void Invert_Blob(Object)
OBJECT *Object;
  {
  ((BLOB *) Object)->Inverted = 1 - ((BLOB *)Object)->Inverted;
  }

void *Copy_Blob (Object)
OBJECT *Object;
  {
  BLOB *New;
  unsigned i, cnt;

  New = Create_Blob();
  *New = * ((BLOB *)Object);

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

  cnt = ((BLOB *)Object)->count;
  New->list = (Blob_Element **)malloc(cnt * sizeof(Blob_Element *));
  if (New->list == NULL) 
    MAError("blob data");
  for (i=0;i<cnt;i++) 
    {
    New->list[i] = (Blob_Element *)malloc(sizeof(Blob_Element));
    if (New->list[i] == NULL)
      MAError("blob data");
    memcpy(New->list[i], ((BLOB *)Object)->list[i], sizeof(Blob_Element));
    }

  New->intervals = (Blob_Interval *)malloc(2*New->count*sizeof(Blob_Interval));
  if (New->intervals == NULL)
    MAError("blob data");

  return (New);
  }

/* Allocate a blob. */
BLOB *Create_Blob()
  {
  BLOB *New;

  if ((New = (BLOB *) malloc (sizeof (BLOB))) == NULL)
    MAError ("blob");

  INIT_OBJECT_FIELDS(New, BLOB_OBJECT, &Blob_Methods)

    New->Trans = NULL;
  New->Inverted = FALSE;
  New->count = 0;
  New->threshold = 0.0;
  New->list = NULL;
  New->intervals = NULL;
  New->Sturm_Flag = FALSE;

  return (New);
  }

void Destroy_Blob (Object)
OBJECT *Object;
  {
  unsigned i;

  Destroy_Transform(((BLOB *)Object)->Trans);
  for (i=0;i < (unsigned)((BLOB *)Object)->count;i++)
    free (((BLOB *)Object)->list[i]);
  free (((BLOB *)Object)->list);
  free (((BLOB *)Object)->intervals);
  free (Object);
  }