/****************************************************************************
 *                   fur.cpp
 *
 *  This module contains all functions for participating fur.
 *
 *  Thomas Willhalm 1999
 *
 * from Persistence of Vision(tm) Ray Tracer version 3.6.
 * Copyright 1991-2003 Persistence of Vision Team
 * Copyright 2003-2004 Persistence of Vision Raytracer Pty. Ltd.
 *---------------------------------------------------------------------------
 * 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 developers. There are strict rules
 * regarding how you are permitted to use this file. These rules are contained
 * in the distribution and derivative versions licenses which should have been
 * provided with this file.
 *
 * These licences may be found online, linked from the end-user license
 * agreement that is located at http://www.povray.org/povlegal.html
 *---------------------------------------------------------------------------
 * 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.
 *---------------------------------------------------------------------------
 *
 *===========================================================================
 * This file is part of MegaPOV, a modified and unofficial version of POV-Ray
 * For more information on MegaPOV visit our website:
 * http://megapov.inetart.net/
 *===========================================================================
 *
 * $RCSfile: fur.cpp,v $
 * $Revision: 1.8 $
 * $Author: chris $
 *
 *****************************************************************************/

/****************************************************************************
* The fur in this picture has been created by combining 
* 
* [1] J.T.Kajiya, T.L.Kay: "Rendering Fur with Three Dimensional Textures" in
* Proc. SIGGRAPH '89, Computer Graphics 23/3, July 1989.
* 
* and
* 
* [2] K.Perlin, E.M.Hoffert: "Hypertexture" in
* Proc. SIGGRAPH '89, Computer Graphics 23/3, July 1989.
* 
* More precisely, the lighting model of [1] was combined with the hair
* distribution of [2]. Some noise and gravity has been added to create
* waves and to let the hairs hang down. The fur actually consists of two
* media containers ("hypertexures") one with thick, short hairs and one with
* longer hair. The container of the latter is modified by a noise function 
* (and therefore is described by an implicit function).
* 
* If you plan to invest some further work in this stuff, please feel free to
* contact me.
* 
* Thomas(at)Willhalm.de
*****************************************************************************/

#include "frame.h"
#include "vector.h"
#include "matrices.h"
#include "texture.h"
#ifdef FUR_PATCH
#include "fur.h"
BEGIN_POV_NAMESPACE

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

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

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

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

static void Project(FUR *Fur, VECTOR TPoint, VECTOR FPoint, VECTOR Tangent);
static void Pertubate(FUR *Fur, VECTOR TPoint, VECTOR FPoint, VECTOR TPoint2);
static DBL bias(DBL t, DBL b);
static DBL gain(DBL t, DBL g);

/*****************************************************************************
*
* FUNCTION
*
*   Create_Fur
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
*   pointer to the created fur
*
* AUTHOR
*
*   Thomas Willhalm (ThW)
*
* DESCRIPTION   : Allocate memory for new fur and initialize it to
*                 system default values.
*
* CHANGES
*
******************************************************************************/

FUR *Create_Fur(void)
{
  FUR* New;
  New  = (FUR *)POV_MALLOC(sizeof (FUR), "fur");
  New->Type                  = SPHERE_FUR;
  New->Vals.Sphere.Radius    = 1.0;
  New->Vals.Sphere.Center[0] = 0.0;
  New->Vals.Sphere.Center[1] = 0.0;
  New->Vals.Sphere.Center[2] = 0.0;
  New->Ratio                 = 0.3;
  New->Falloff               = 0.9;
  New->Frequency             = 200.0;
  New->Force[0]              = 0.0;
  New->Force[1]              = 0.0;
  New->Force[2]              = 0.0;
  New->Diffuse               = 5.0;
  New->Reflection            = 1.0;
  New->Reflection_Exponent   = 1.0;
  New->Waves_Frequency       = 0.0;
  New->Waves_Factor          = 0.0;
  New->Trans                 = NULL;
  return New;
}

/*****************************************************************************
*
* FUNCTION
*
*   Copy_Fur
*
* INPUT
*
*   Old -- point to fur to be copied
*
* RETURNS
*
*   pointer to the created fur
*
* AUTHOR
*
*   Thomas Willhalm (ThW)
*
* DESCRIPTION   : Allocate memory for new pattern and initialize it to
*                 values in existing pattern Old.
*
* CHANGES
*
******************************************************************************/

FUR *Copy_Fur(FUR *Old)
{
  FUR* New;

  if (Old != NULL)
  {
    New = Create_Fur();
    *New = *Old;
    New->Trans = Copy_Transform(Old->Trans);
  }
  else
  {
    New = NULL;
  }

  return New;
}

/*****************************************************************************
*
* FUNCTION
*
*   Destroy_Fur
*
* INPUT
*
*   Fur - pointer to fur to destroy
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   Thomas Willhalm (ThW)
*
* DESCRIPTION   : free all memory associated with given fur
*
* CHANGES
*
******************************************************************************/

void Destroy_Fur(FUR *Fur)
{
  if (Fur != NULL)
  {
    Destroy_Transform(Fur->Trans);
    POV_FREE(Fur);
  }
}

/*****************************************************************************
*
* FUNCTION
*
*   Transform_Fur
*
* INPUT
*
*   Fur - pointer to fur to transform
*   Trans - pointer to transformation struct
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   Thomas Willhalm (ThW)
*
* DESCRIPTION   :  Transform fur
*
* CHANGES
*
******************************************************************************/

void Transform_Fur(FUR *Fur, TRANSFORM *Trans)
{
  if(Fur!=NULL)
  {
    if (Fur->Trans==NULL)
    {
      Fur->Trans = Create_Transform();
    }

    Compose_Transforms(Fur->Trans, Trans);
    MTransDirection(Fur->Force, Fur->Force, Trans);
  }
}

/*****************************************************************************
*
* FUNCTION
*
*   evaluate_fur_density
*
* INPUT
*
*  Fur -  pointer to fur struct
*  EPoint - location to evaluate
*
* OUTPUT
*
*  Tangent - tangent at location EPoint
*
* RETURNS
*
*   density value
*
* AUTHOR
*
*   Thomas Willhalm (ThW)
*
* DESCRIPTION   :
*
*   Calculates the density at EPoint as well as a vector that is the direction
*   of a tangent of the hairs.
*
* CHANGES
*
******************************************************************************/

DBL evaluate_fur_density(FUR *Fur, VECTOR EPoint, VECTOR Tangent)
{
  DBL s;
  VECTOR TPoint;  /* transformed */
  VECTOR FPoint;  /* transformed and projected */
  VECTOR TPointPertubated; /* transformed and pertubated */
  VECTOR FPointPertubated; /* transformed, pertubated and projected */
  VECTOR TPoint2; /* slightly moved TPoint for calculation of Tangent */
  VECTOR FPoint2;
  VECTOR TPoint2Pertubated;
  VECTOR Tangent2; /* Tangent without pertubation */

  /* transform EPoint */

  if (Fur->Trans==NULL)
    { Assign_Vector(TPoint, EPoint); }
  else
    { MInvTransPoint(TPoint, EPoint, Fur->Trans); }

  Project(Fur, TPoint, FPoint, Tangent2);
  Pertubate(Fur, TPoint, FPoint, TPointPertubated);

  /* calculate density */
  Project(Fur, TPointPertubated, FPointPertubated, Tangent); /* Tangent is dummy */
  VScaleEq(FPointPertubated, Fur->Frequency);
  s = gain(bias(Noise(FPointPertubated,NULL), Fur->Ratio), Fur->Falloff);

  /* calculate the tangent */
  VAddScaled(TPoint2, TPoint, 0.1, Tangent2);
  Project(Fur, TPoint2, FPoint2, Tangent2); /* Tangent2 is dummy */
  Pertubate(Fur, TPoint2, FPoint2, TPoint2Pertubated);
  VSub(Tangent, TPoint2Pertubated, TPointPertubated);
  VNormalizeEq(Tangent);

  return s;
}

/*****************************************************************************
*
* FUNCTION
*
*   Pertubate
*
* INPUT
*
*  TPoint  - vector to pertubate
*  FPoint  - projection of TPoint on surface
*
* OUTPUT
*
*  TPoint2 - pertubated vector
*
* RETURNS
*
* AUTHOR
*
*   Thomas Willhalm (ThW)
*
* DESCRIPTION   :
*
*   Applies force and waves.
*
* CHANGES
*
******************************************************************************/


static void Pertubate(FUR *Fur, VECTOR TPoint, VECTOR FPoint, VECTOR TPoint2)
{
  VECTOR FPoint2;
  DBL d;

  /* apply force */
  VDist(d, TPoint, FPoint);
  VAddScaled(TPoint2, TPoint, (-d*d), Fur->Force);

  /* apply waves */
  VScale(FPoint2,FPoint, Fur->Waves_Frequency);
  DNoise(FPoint2,FPoint2);
  VAddScaledEq(TPoint2, d*d*(Fur->Waves_Factor), FPoint2);
}

/*****************************************************************************
*
* FUNCTION
*
*   Project
*
* INPUT
*
*  Fur -  pointer to fur struct
*  TPoint - location to evaluate (already transformed)
*
* OUTPUT
*
*  FPoint  - projection of TPoint on surface
*  Tangent - tangent vector of fur (orthogonal to surface)
*
* RETURNS
*
* AUTHOR
*
*   Thomas Willhalm (ThW)
*
* DESCRIPTION   :
*
*   Projects TPoint on the surface of the underlaying structure.
*
* CHANGES
*
******************************************************************************/

static void Project_Smooth_Triangle(FUR *Fur, VECTOR TPoint, VECTOR FPoint, VECTOR Tangent)
{
  int Axis;
  DBL u, v, t, s;
  VECTOR PP1, PP2, PP3, TPointMinusPP1;

  VDot(t, TPoint, Fur->Vals.Smooth_Triangle.Normal_Vector);
  t -= Fur->Vals.Smooth_Triangle.Distance;

  VAddScaled(PP1, Fur->Vals.Smooth_Triangle.P1, t, Fur->Vals.Smooth_Triangle.N1);
  VAddScaled(PP2, Fur->Vals.Smooth_Triangle.P2, t, Fur->Vals.Smooth_Triangle.N2);
  VAddScaled(PP3, Fur->Vals.Smooth_Triangle.P3, t, Fur->Vals.Smooth_Triangle.N3);

  VSub(TPointMinusPP1, TPoint, PP1);

  VDot(u, TPointMinusPP1, Fur->Vals.Smooth_Triangle.Perp); /* not accurate */

  if (u < EPSILON)
  {
    Assign_Vector(Tangent, Fur->Vals.Smooth_Triangle.N1);
  }
  else
  {

    Axis = Fur->Vals.Smooth_Triangle.vAxis;

    v = (TPointMinusPP1[Axis]/u + PP1[Axis]  - PP2[Axis]) / (PP3[Axis]-PP2[Axis]);


    Tangent[X] = Fur->Vals.Smooth_Triangle.N1[X]
                 + u*(Fur->Vals.Smooth_Triangle.N2[X]-Fur->Vals.Smooth_Triangle.N1[X]
		 + v*(Fur->Vals.Smooth_Triangle.N3[X]-Fur->Vals.Smooth_Triangle.N2[X]));
    Tangent[Y] = Fur->Vals.Smooth_Triangle.N1[Y]
                 + u*(Fur->Vals.Smooth_Triangle.N2[Y]-Fur->Vals.Smooth_Triangle.N1[Y]
		 + v*(Fur->Vals.Smooth_Triangle.N3[Y]-Fur->Vals.Smooth_Triangle.N2[Y]));
    Tangent[Z] = Fur->Vals.Smooth_Triangle.N1[Z]
                 + u*(Fur->Vals.Smooth_Triangle.N2[Z]-Fur->Vals.Smooth_Triangle.N1[Z]
		 + v*(Fur->Vals.Smooth_Triangle.N3[Z]-Fur->Vals.Smooth_Triangle.N2[Z]));
  }
  VNormalizeEq(Tangent);

  VDot(s, TPoint, Fur->Vals.Smooth_Triangle.Normal_Vector);
  s -= Fur->Vals.Smooth_Triangle.Distance;
  VDot(t, Tangent, Fur->Vals.Smooth_Triangle.Normal_Vector);
  s /= -t;
  VAddScaled(FPoint, TPoint, s, Tangent);
}

static void Project(FUR *Fur, VECTOR TPoint, VECTOR FPoint, VECTOR Tangent)
{
  register DBL t;
  switch (Fur->Type)
    {
    case (SPHERE_FUR):
      VSub(Tangent, TPoint, Fur->Vals.Sphere.Center);
      VNormalizeEq(Tangent);
      VAddScaled(FPoint, Fur->Vals.Sphere.Center,
		 Fur->Vals.Sphere.Radius, Tangent);
      break;
    case (TORUS_FUR):
      t = Fur->Vals.Torus.Major_Radius/sqrt(TPoint[0]*TPoint[0] + TPoint[2]*TPoint[2]);
      FPoint[0] = TPoint[0]*t;
      FPoint[1] = 0;
      FPoint[2] = TPoint[2]*t;
      VSub(Tangent, TPoint, FPoint);
      VNormalizeEq(Tangent);
      VAddScaledEq(FPoint, Fur->Vals.Torus.Minor_Radius, Tangent);
      break;
    case (BOX_FUR):
      VNormalize(Tangent, TPoint);
      VEvaluate(FPoint, Tangent, Fur->Vals.Box.Size);
      break;
    case (PLANE_FUR):
      Tangent[0] = Fur->Vals.Plane.Normal[0];
      Tangent[1] = Fur->Vals.Plane.Normal[1];
      Tangent[2] = Fur->Vals.Plane.Normal[2];
      VDot(t, TPoint, Fur->Vals.Plane.Normal);
      t -= Fur->Vals.Plane.Distance;
      VSubScaled(FPoint, TPoint, t, Fur->Vals.Plane.Normal);
      break;
    case (SMOOTH_TRIANGLE_FUR):
      Project_Smooth_Triangle(Fur, TPoint, FPoint, Tangent);
      break;
    case (CYLINDER_FUR):
      VSub(Tangent, TPoint, Fur->Vals.Cylinder.Point);
      VDot(t, Tangent, Fur->Vals.Cylinder.Direction);
      if (t>0) { /* beyond the front of the cylinder */
	VNormalizeEq(Tangent);
	VScale(FPoint, Tangent, Fur->Vals.Cylinder.Radius);
	VAddEq(FPoint, Fur->Vals.Cylinder.Point);
      } else {
	if (t>Fur->Vals.Cylinder.Length) /* beyond the back of the cylinder */
	  { t = Fur->Vals.Cylinder.Length; }

	VSubScaledEq(Tangent, t, Fur->Vals.Cylinder.Direction);
	VNormalizeEq(Tangent);
	VScale(FPoint, Tangent, Fur->Vals.Cylinder.Radius);
	VAddScaledEq(FPoint, t, Fur->Vals.Cylinder.Direction);
	VAddEq(FPoint, Fur->Vals.Cylinder.Point);
      }
	break;
    }
}

static DBL bias(DBL t, DBL b)
{
   return pow(t,log(b)/log(0.5));
}

static DBL gain(DBL t, DBL g)
{
   if (t<0.5)
        return bias(2*t,1-g)/2;
   else
        return 1-bias(2-2*t,1-g)/2;
}

END_POV_NAMESPACE

#endif
