/**************************************************************************
*                pattern.c
*
*  This module implements texturing functions that return a value to be
*  used in a pigment or normal.
*
*  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.
*
*****************************************************************************/

/*
 * Some texture ideas garnered from SIGGRAPH '85 Volume 19 Number 3,
 * "An Image Synthesizer" By Ken Perlin.
 * Further Ideas Garnered from "The RenderMan Companion" (Addison Wesley).
 */

#include "frame.h"
#include "vector.h"
#include "povproto.h"
#include "matrices.h"
#include "pattern.h"
#include "texture.h"
#include "image.h"
#include "txttest.h"
#include "colour.h"



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


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

static DBL agate PARAMS((VECTOR EPoint, TPATTERN *TPat));
static DBL brick PARAMS((VECTOR EPoint, TPATTERN *TPat));
static DBL checker PARAMS((VECTOR EPoint));
static DBL crackle PARAMS((VECTOR EPoint));
static DBL gradient PARAMS((VECTOR EPoint, TPATTERN *TPat));
static DBL granite PARAMS((VECTOR EPoint));
static DBL leopard PARAMS((VECTOR EPoint));
static DBL mandel PARAMS((VECTOR EPoint, TPATTERN *TPat));
static DBL marble PARAMS((VECTOR EPoint, TPATTERN *TPat));
static DBL onion PARAMS((VECTOR EPoint));
static DBL radial PARAMS((VECTOR EPoint));
static DBL spiral1 PARAMS((VECTOR EPoint, TPATTERN *TPat));
static DBL spiral2 PARAMS((VECTOR EPoint, TPATTERN *TPat));
static DBL wood PARAMS((VECTOR EPoint, TPATTERN *TPat));
static DBL hexagon PARAMS((VECTOR EPoint));
static long PickInCube PARAMS((VECTOR tv, VECTOR p1));

static DBL ripples_pigm PARAMS((VECTOR EPoint, TPATTERN *TPat));
static DBL waves_pigm  PARAMS((VECTOR EPoint, TPATTERN *TPat));
static DBL dents_pigm  PARAMS((VECTOR EPoint));
static DBL wrinkles_pigm PARAMS((VECTOR EPoint));
static DBL quilted_pigm PARAMS((VECTOR EPoint, TPATTERN *TPat));
static TURB *Search_For_Turb PARAMS((WARP *Warps));
/* static TURB *Copy_Turb PARAMS((TURB *Old));   Unused function [AED] */



/*****************************************************************************
*
* FUNCTION
*
*   agate
*
* INPUT
*
*   EPoint -- The point in 3d space at which the pattern is evaluated.
*   TPat   -- Texture pattern struct
*
* OUTPUT
*
* RETURNS
*
*   DBL value in the range 0.0 to 1.0
*
* AUTHOR
*
*   POV-Ray Team
*
* DESCRIPTION
*
* CHANGES
*   Oct 1994    : adapted from agate pigment by [CY]
*
******************************************************************************/

static DBL agate (EPoint, TPat)
VECTOR EPoint;
TPATTERN *TPat;
{
  register DBL noise, turb_val;
  TURB* Turb;

  Turb=Search_For_Turb(TPat->Warps);

  turb_val = TPat->Vals.Agate_Turb_Scale * Turbulence(EPoint,Turb);

  noise = 0.5 * (cycloidal(1.3 * turb_val + 1.1 * EPoint[Z]) + 1.0);

  if (noise < 0.0)
  {
    noise = 0.0;
  }
  else
  {
    noise = min(1.0, noise);
    noise = pow(noise, 0.77);
  }

  return(noise);
}


/*****************************************************************************
*
* FUNCTION
*
*   brick
*
* INPUT
*
*   EPoint -- The point in 3d space at which the pattern
*   is evaluated.
*   TPat   -- Texture pattern struct
*   
* OUTPUT
*   
* RETURNS
*
*   DBL value exactly 0.0 or 1.0
*   
* AUTHOR
*
*   Dan Farmer
*   
* DESCRIPTION
*
* CHANGES
*   Oct 1994    : adapted from pigment by [CY]
*
******************************************************************************/

static DBL brick (EPoint, TPat)
VECTOR EPoint;
TPATTERN *TPat;
{
  int ibrickx, ibricky, ibrickz;
  DBL brickheight, brickwidth, brickdepth;
  DBL brickmortar, mortarheight, mortarwidth, mortardepth;
  DBL brickx, bricky, brickz;
  DBL x, y, z, fudgit;

  fudgit=Small_Tolerance+TPat->Vals.Brick.Mortar;

  x =  EPoint[X]+fudgit;
  y =  EPoint[Y]+fudgit;
  z =  EPoint[Z]+fudgit;

  brickwidth  = TPat->Vals.Brick.Size[X];
  brickheight = TPat->Vals.Brick.Size[Y];
  brickdepth  = TPat->Vals.Brick.Size[Z];
  brickmortar = (DBL)TPat->Vals.Brick.Mortar;

  mortarwidth  = brickmortar / brickwidth;
  mortarheight = brickmortar / brickheight;
  mortardepth  = brickmortar / brickdepth;

  /* 1) Check mortar layers in the X-Z plane (ie: top view) */

  bricky = y / brickheight;
  ibricky = (int) bricky;
  bricky -= (DBL) ibricky;

  if (bricky < 0.0)
  {
    bricky += 1.0;
  }

  if (bricky <= mortarheight)
  {
    return(0.0);
  }

  bricky = (y / brickheight) * 0.5;
  ibricky = (int) bricky;
  bricky -= (DBL) ibricky;

  if (bricky < 0.0)
  {
    bricky += 1.0;
  }


  /* 2) Check ODD mortar layers in the Y-Z plane (ends) */

  brickx = (x / brickwidth);
  ibrickx = (int) brickx;
  brickx -= (DBL) ibrickx;

  if (brickx < 0.0)
  {
    brickx += 1.0;
  }

  if ((brickx <= mortarwidth) && (bricky <= 0.5))
  {
    return(0.0);
  }

  /* 3) Check EVEN mortar layers in the Y-Z plane (ends) */

  brickx = (x / brickwidth) + 0.5;
  ibrickx = (int) brickx;
  brickx -= (DBL) ibrickx;

  if (brickx < 0.0)
  {
    brickx += 1.0;
  }

  if ((brickx <= mortarwidth) && (bricky > 0.5))
  {
    return(0.0);
  }

  /* 4) Check ODD mortar layers in the Y-X plane (facing) */

  brickz = (z / brickdepth);
  ibrickz = (int) brickz;
  brickz -= (DBL) ibrickz;

  if (brickz < 0.0)
  {
    brickz += 1.0;
  }

  if ((brickz <= mortardepth) && (bricky > 0.5))
  {
    return(0.0);
  }

  /* 5) Check EVEN mortar layers in the X-Y plane (facing) */

  brickz = (z / brickdepth) + 0.5;
  ibrickz = (int) brickz;
  brickz -= (DBL) ibrickz;

  if (brickz < 0.0)
  {
    brickz += 1.0;
  }

  if ((brickz <= mortardepth) && (bricky <= 0.5))
  {
    return(0.0);
  }

  /* If we've gotten this far, color me brick. */

  return(1.0);
}


/*****************************************************************************
*
* FUNCTION
*
*   checker
*
* INPUT
*
*   EPoint -- The point in 3d space at which the pattern
*   is evaluated.
*
* OUTPUT
*
* RETURNS
*
*   DBL value exactly 0.0 or 1.0
*
* AUTHOR
*
*   POV-Team
*
* DESCRIPTION
*
* CHANGES
*   Oct 1994    : adapted from pigment by [CY]
*
******************************************************************************/

static DBL checker (EPoint)
VECTOR EPoint;
{
  int value;

  value = (int)(floor(EPoint[X]+Small_Tolerance) +
                floor(EPoint[Y]+Small_Tolerance) +
                floor(EPoint[Z]+Small_Tolerance));

  if (value & 1)
  {
    return (1.0);
  }
  else
  {
    return (0.0);
  }
}



/*****************************************************************************
*
* FUNCTION
*
*   crackle
*
* INPUT
*
*   EPoint -- The point in 3d space at which the pattern
*   is evaluated.
* OUTPUT
*
* RETURNS
*
*   DBL value in the range 0.0 to 1.0
*
* AUTHOR
*
*   Jim McElhiney
*
* DESCRIPTION
*
*   "crackle":
*
*   New colour function by Jim McElhiney,
*     CompuServe 71201,1326, aka mcelhiney@acm.org
*
*   Large scale, without turbulence, makes a pretty good stone wall.
*   Small scale, without turbulence, makes a pretty good crackle ceramic glaze.
*   Highly turbulent (with moderate displacement) makes a good marble, solving
*   the problem of apparent parallel layers in Perlin's method.
*   2 octaves of full-displacement turbulence make a great "drizzled paint"
*   pattern, like a 1950's counter top.
*   Rule of thumb:  put a single colour transition near 0 in your colour map.
*
*   Mathematically, the set crackle(p)=0 is a 3D Voronoi diagram of a field of
*   semirandom points, and crackle(p)>0 is distance from set along shortest path.
*   (A Voronoi diagram is the locus of points equidistant from their 2 nearest
*   neighbours from a set of disjoint points, like the membranes in suds are
*   to the centres of the bubbles).
*
*   All "crackle" specific source code and examples are in the public domain.
*
* CHANGES
*   Oct 1994    : adapted from pigment by [CY]
*
******************************************************************************/

static DBL crackle (EPoint)
VECTOR EPoint;
{
  int    i;
  long   thisseed;
  DBL    sum, minsum, minsum2, tf;
  VECTOR sv, tv, dv, t1, add;

  static int cvc;
  static long lastseed = 0x80000000;
  static VECTOR cv[81];

  Assign_Vector(tv,EPoint);

  /*
   * Check to see if the input point is in the same unit cube as the last
   * call to this function, to use cache of cubelets for speed.
   */

  thisseed = PickInCube(tv, t1);

  if (thisseed != lastseed)
  {
    /*
     * No, not same unit cube.  Calculate the random points for this new
     * cube and its 80 neighbours which differ in any axis by 1 or 2.
     * Why distance of 2?  If there is 1 point in each cube, located
     * randomly, it is possible for the closest random point to be in the
     * cube 2 over, or the one two over and one up.  It is NOT possible
     * for it to be two over and two up.  Picture a 3x3x3 cube with 9 more
     * cubes glued onto each face.
     */

    /* Now store a points for this cube and each of the 80 neighbour cubes. */

    cvc = 0;

    for (add[X] = -2.0; add[X] < 2.5; add[X] +=1.0)
    {
      for (add[Y] = -2.0; add[Y] < 2.5; add[Y] += 1.0)
      {
        for (add[Z] = -2.0; add[Z] < 2.5; add[Z] += 1.0)
        {
          /* For each cubelet in a 5x5 cube. */

          if ((fabs(add[X])>1.5)+(fabs(add[Y])>1.5)+(fabs(add[Z])>1.5) <= 1.0)
          {
            /* Yes, it's within a 3d knight move away. */

            VAdd(sv, tv, add);

            PickInCube(sv, t1);

            cv[cvc][X] = t1[X];
            cv[cvc][Y] = t1[Y];
            cv[cvc][Z] = t1[Z];
            cvc++;
          }
        }
      }
    }

    lastseed = thisseed;
  }

  /*
   * Find the 2 points with the 2 shortest distances from the input point.
   * Loop invariant:  minsum is shortest dist, minsum2 is 2nd shortest
   */

  /* Set up the loop so the invariant is true:  minsum <= minsum2 */

  VSub(dv, cv[0], tv);  minsum  = VSumSqr(dv);
  VSub(dv, cv[1], tv);  minsum2 = VSumSqr(dv);

  if (minsum2 < minsum)
  {
    tf = minsum; minsum = minsum2; minsum2 = tf;
  }

  /* Loop for the 81 cubelets to find closest and 2nd closest. */

  for (i = 2; i < cvc; i++)
  {
    VSub(dv, cv[i], tv);

    sum = VSumSqr(dv);

    if (sum < minsum)
    {
      minsum2 = minsum;
      minsum = sum;
    }
    else
    {
      if (sum < minsum2)
      {
        minsum2 = sum;
      }
    }
  }

  /* Crackle value is absolute value of diff in dist to closest 2 points. */

  tf = sqrt(minsum2) - sqrt(minsum);      /* minsum is known <= minsum2 */

  /*
   * Note that the theoretical range of this function is 0 to root 3.
   * In practice, it rarely exceeds 0.9, and only very rarely 1.0
   */

  return min(tf, 1.);
}


/*****************************************************************************
*
* FUNCTION
*
*   gradient
*
* INPUT
*
*   EPoint -- The point in 3d space at which the pattern
*   is evaluated.
*   
* OUTPUT
*   
* RETURNS
*
*   DBL value in the range 0.0 to 1.0
*   
* AUTHOR
*
*   POV-Ray Team
*   
* DESCRIPTION
*
*   Gradient Pattern - gradient based on the fractional values of
*   x, y or z, based on whether or not the given directional vector is
*   a 1.0 or a 0.0.
*   The basic concept of this is from DBW Render, but Dave Wecker's
*   only supports simple Y axis gradients.
*
* CHANGES
*   Oct 1994    : adapted from pigment by [CY]
*
******************************************************************************/

static DBL gradient (EPoint, TPat)
VECTOR EPoint;
TPATTERN *TPat;
{
  register int i;
  register DBL temp;
  DBL value = 0.0;

  for (i=X; i<=Z; i++)
  {
    if (TPat->Vals.Gradient[i] != 0.0)
    {
      temp = fabs(EPoint[i]);

      value += fmod(temp,1.0);
    }
  }

  /* Clamp to 1.0. */

  value = ((value > 1.0) ? fmod(value, 1.0) : value);

  return(value);
}



/*****************************************************************************
*
* FUNCTION
*
*   granite
*
* INPUT
*
*   EPoint -- The point in 3d space at which the pattern
*   is evaluated.
*   
* OUTPUT
*   
* RETURNS
*
*   DBL value in the range 0.0 to 1.0
*   
* AUTHOR
*
*   POV-Ray Team
*   
* DESCRIPTION
*
*   Granite - kind of a union of the "spotted" and the "dented" textures,
*   using a 1/f fractal noise function for color values. Typically used
*   with small scaling values. Should work with colour maps for pink granite.
*
* CHANGES
*   Oct 1994    : adapted from pigment by [CY]
*
******************************************************************************/

static DBL granite (EPoint)
VECTOR EPoint;
{
  register int i;
  register DBL temp, noise = 0.0, freq = 1.0;
  VECTOR tv1,tv2;

  VScale(tv1,EPoint,4.0);

  for (i = 0; i < 6 ; freq *= 2.0, i++)
  {
    VScale(tv2,tv1,freq);
    temp = 0.5 - Noise (tv2);

    temp = fabs(temp);

    noise += temp / freq;
  }

  return(noise);
}



/*****************************************************************************
*
* FUNCTION
*
*   leopard
*
* INPUT
*
*   EPoint -- The point in 3d space at which the pattern
*   is evaluated.
*
* OUTPUT
*
* RETURNS
*
*   DBL value in the range 0.0 to 1.0
*
* AUTHOR
*
*   Scott Taylor
*
* DESCRIPTION
*
* CHANGES
*   Jul 1991 : Creation.
*   Oct 1994 : adapted from pigment by [CY]
*
******************************************************************************/

static DBL leopard (EPoint)
VECTOR EPoint;
{
  register DBL value, temp1, temp2, temp3;

  /* This form didn't work with Zortech 386 compiler */
  /* value = Sqr((sin(x)+sin(y)+sin(z))/3); */
  /* So we break it down. */

  temp1 = sin(EPoint[X]);
  temp2 = sin(EPoint[Y]);
  temp3 = sin(EPoint[Z]);

  value = Sqr((temp1 + temp2 + temp3) / 3.0);

  return(value);
}



/*****************************************************************************
*
* FUNCTION
*
*   mandel
*
* INPUT
*
*   EPoint -- The point in 3d space at which the pattern
*   is evaluated.
*
* OUTPUT
*
* RETURNS
*
*   DBL value in the range 0.0 to 1.0
*
* AUTHOR
*
*   submitted by user, name lost (sorry)
*
* DESCRIPTION
* The mandel pattern computes the standard Mandelbrot fractal pattern and
* projects it onto the X-Y plane.  It uses the X and Y coordinates to compute
* the Mandelbrot set.
*
* CHANGES
*   Oct 1994 : adapted from pigment by [CY]
*
******************************************************************************/

static DBL mandel (EPoint, TPat)
VECTOR EPoint;
TPATTERN *TPat;
{
  int it_max, col;
  DBL a, b, cf, a2, b2, x, y;

  a = x = EPoint[X]; a2 = Sqr(a);
  b = y = EPoint[Y]; b2 = Sqr(b);

  it_max = TPat->Vals.Iterations;

  for (col = 0; col < it_max; col++)
  {
    b  = 2.0 * a * b + y;
    a  = a2 - b2 + x;

    a2 = Sqr(a);
    b2 = Sqr(b);

    if ((a2 + b2) > 4.0)
    {
      break;
    }
  }

  cf = (DBL)col / (DBL)it_max;

  return(cf);
}



/*****************************************************************************
*
* FUNCTION
*
*   marble
*
* INPUT
*
*   EPoint -- The point in 3d space at which the pattern
*   is evaluated.
*   TPat   -- Texture pattern struct
*
* OUTPUT
*
* RETURNS
*
*   DBL value in the range 0.0 to 1.0
*
* AUTHOR
*
*   POV-Ray Team
*
* DESCRIPTION
*
* CHANGES
*   Oct 1994 : adapted from pigment by [CY]
*
******************************************************************************/

static DBL marble (EPoint, TPat)
VECTOR EPoint;
TPATTERN *TPat;
{
  register DBL turb_val;
  TURB *Turb;

  if ((Turb=Search_For_Turb(TPat->Warps)) != NULL)
  {
    turb_val = Turb->Turbulence[X] * Turbulence(EPoint,Turb);
  }
  else
  {
    turb_val = 0.0;
  }

  return(EPoint[X] + turb_val);
}



/*****************************************************************************
*
* FUNCTION
*
*   onion
*
* INPUT
*
*   EPoint -- The point in 3d space at which the pattern
*   is evaluated.
*
* OUTPUT
*
* RETURNS
*
*   DBL value in the range 0.0 to 1.0
*
* AUTHOR
*
*   Scott Taylor
*
* DESCRIPTION
*
* CHANGES
*   Jul 1991 : Creation.
*   Oct 1994 : adapted from pigment by [CY]
*
******************************************************************************/

static DBL onion (EPoint)
VECTOR EPoint;
{
  /* The variable noise is not used as noise in this function */

  register DBL noise;

/*
   This ramp goes 0-1,1-0,0-1,1-0...

   noise = (fmod(sqrt(Sqr(x)+Sqr(y)+Sqr(z)),2.0)-1.0);

   if (noise<0.0) {noise = 0.0-noise;}
*/

  /* This ramp goes 0-1, 0-1, 0-1, 0-1 ... */

  noise = (fmod(sqrt(Sqr(EPoint[X])+Sqr(EPoint[Y])+Sqr(EPoint[Z])), 1.0));

  return(noise);
}



/*****************************************************************************
*
* FUNCTION
*
*   radial
*
* INPUT
*
*   EPoint -- The point in 3d space at which the pattern
*   is evaluated.
*
* OUTPUT
*
* RETURNS
*
*   DBL value in the range 0.0 to 1.0
*
* AUTHOR
*
*   Chris Young -- new in vers 2.0
*
* DESCRIPTION
*
* CHANGES
*   Oct 1994 : adapted from pigment by [CY]
*
******************************************************************************/

static DBL radial (EPoint)
VECTOR EPoint;
{
  register DBL value;

  if ((fabs(EPoint[X])<0.001) && (fabs(EPoint[Z])<0.001))
  {
    value = 0.25;
  }
  else
  {
    value = 0.25 + (atan2(EPoint[X],EPoint[Z]) + M_PI) / TWO_M_PI;
  }

  return(value);
}



/*****************************************************************************
*
* FUNCTION
*
*   spiral1
*
* INPUT
*
*   EPoint -- The point in 3d space at which the pattern
*   is evaluated.
*   TPat   -- Texture pattern struct
*
* OUTPUT
*
* RETURNS
*
*   DBL value in the range 0.0 to 1.0
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*   Spiral whirles around z-axis.
*   The number of "arms" is defined in the TPat.
*
* CHANGES
*   Aug 1994 : Creation.
*   Oct 1994 : adapted from pigment by [CY]
*
******************************************************************************/

static DBL spiral1(EPoint, TPat)
VECTOR EPoint;
TPATTERN *TPat;
{
  DBL rad, phi, turb_val;
  DBL x = EPoint[X];
  DBL y = EPoint[Y];
  DBL z = EPoint[Z];
  TURB *Turb;

  if ((Turb=Search_For_Turb(TPat->Warps)) != NULL)
  {
    turb_val = Turb->Turbulence[X] * Turbulence(EPoint,Turb);
  }
  else
  {
    turb_val = 0.0;
  }

  /* Get distance from z-axis. */

  rad = sqrt(x * x + y * y);

  /* Get angle in x,y-plane (0...2 PI). */

  if (rad == 0.0)
  {
    phi = 0.0;
  }
  else
  {
    if (x < 0.0)
    {
      phi = 3.0 * M_PI_2 - asin(y / rad);
    }
    else
    {
      phi = M_PI_2 + asin(y / rad);
    }
  }

  return(z + rad + (DBL)TPat->Vals.Arms * phi / TWO_M_PI + turb_val);
}



/*****************************************************************************
*
* FUNCTION
*
*   spiral2
*
* INPUT
*
*   EPoint -- The point in 3d space at which the pattern
*   is evaluated.
*   TPat   -- Texture pattern struct
*
* OUTPUT
*
* RETURNS
*
*   DBL value in the range 0.0 to 1.0
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*   Spiral whirles around z-axis.
*   The number of "arms" is defined in the TPat.
*
* CHANGES
*   Aug 1994 : Creation.
*   Oct 1994 : adapted from pigment by [CY]
*
******************************************************************************/


static DBL spiral2(EPoint, TPat)
VECTOR EPoint;
TPATTERN *TPat;
{
  DBL rad, phi, turb_val;
  DBL x = EPoint[X];
  DBL y = EPoint[Y];
  DBL z = EPoint[Z];
  TURB *Turb;

  if ((Turb=Search_For_Turb(TPat->Warps)) != NULL)
  {
    turb_val = Turb->Turbulence[X] * Turbulence(EPoint,Turb);
  }
  else
  {
    turb_val = 0.0;
  }

  /* Get distance from z-axis. */

  rad = sqrt(x * x + y * y);

  /* Get angle in x,y-plane (0...2 PI) */

  if (rad == 0.0)
  {
    phi = 0.0;
  }
  else
  {
    if (x < 0.0)
    {
      phi = 3.0 * M_PI_2 - asin(y / rad);
    }
    else
    {
      phi = M_PI_2 + asin(y / rad);
    }
  }

  turb_val = Triangle_Wave(z + rad + (DBL)TPat->Vals.Arms * phi / TWO_M_PI +
                           turb_val);

  return(Triangle_Wave(rad) + turb_val);
}



/*****************************************************************************
*
* FUNCTION
*
*   wood
*
* INPUT
*
*   EPoint -- The point in 3d space at which the pattern
*   is evaluated.
*   TPat   -- Texture pattern struct
*
* OUTPUT
*
* RETURNS
*
*   DBL value in the range 0.0 to 1.0
*
* AUTHOR
*
*   POV-Ray Team
*
* DESCRIPTION
*
* CHANGES
*   Oct 1994 : adapted from pigment by [CY]
*
******************************************************************************/


static DBL wood (EPoint, TPat)
VECTOR EPoint;
TPATTERN *TPat;
{
  register DBL length;
  VECTOR WoodTurbulence;
  VECTOR point;
  DBL x=EPoint[X];
  DBL y=EPoint[Y];
  TURB *Turb;

  if ((Turb=Search_For_Turb(TPat->Warps)) != NULL)
  {
    DTurbulence (WoodTurbulence, EPoint,Turb);
    point[X] = cycloidal((x + WoodTurbulence[X]) * Turb->Turbulence[X]);
    point[Y] = cycloidal((y + WoodTurbulence[Y]) * Turb->Turbulence[Y]);
  }
  else
  {
    point[X] = 0.0;
    point[Y] = 0.0;
  }
  point[Z] = 0.0;

  point[X] += x;
  point[Y] += y;

  /* point[Z] += z; Deleted per David Buck --  BP 7/91 */

  VLength (length, point);

  return(length);
}


/*****************************************************************************
*
* FUNCTION
*
*   hexagon
*
* INPUT
*
*   EPoint -- The point in 3d space at which the pattern
*   is evaluated.
*
* OUTPUT
*
* RETURNS
*
*   DBL value exactly 0.0, 1.0 or 2.0
*
* AUTHOR
*
*   Ernest MacDougal Campbell III
*   
* DESCRIPTION
*
*   TriHex pattern -- Ernest MacDougal Campbell III (EMC3) 11/23/92
*
*   Creates a hexagon pattern in the XZ plane.
*
*   This algorithm is hard to explain.  First it scales the point to make
*   a few of the later calculations easier, then maps some points to be
*   closer to the Origin.  A small area in the first quadrant is subdivided
*   into a 6 x 6 grid.  The position of the point mapped into that grid
*   determines its color.  For some points, just the grid location is enough,
*   but for others, we have to calculate which half of the block it's in
*   (this is where the atan2() function comes in handy).
*
* CHANGES
*   Nov 1992 : Creation.
*   Oct 1994 : adapted from pigment by [CY]
*
******************************************************************************/

#define xfactor 0.5;         /* each triangle is split in half for the grid */
#define zfactor 0.866025404; /* sqrt(3)/2 -- Height of an equilateral triangle */

static DBL hexagon (EPoint)
VECTOR EPoint;
{
  int xm, zm;
  int brkindx;
  DBL xs, zs, xl, zl, value = 0.0;
  DBL x=EPoint[X];
  DBL z=EPoint[Z];


  /* Keep all numbers positive.  Also, if z is negative, map it in such a
   * way as to avoid mirroring across the x-axis.  The value 5.196152424
   * is (sqrt(3)/2) * 6 (because the grid is 6 blocks high)
   */

  x = fabs(x);

  /* Avoid mirroring across x-axis. */

  z = z < 0.0 ? 5.196152424 - fabs(z) : z;

  /* Scale point to make calcs easier. */

  xs = x/xfactor;
  zs = z/zfactor;

  /* Map points into the 6 x 6 grid where the basic formula works. */

  xs -= floor(xs/6.0) * 6.0;
  zs -= floor(zs/6.0) * 6.0;

  /* Get a block in the 6 x 6 grid. */

  xm = (int) FLOOR(xs) % 6;
  zm = (int) FLOOR(zs) % 6;

  switch (xm)
  {
    /* These are easy cases: Color depends only on xm and zm. */

    case 0:
    case 5:

      switch (zm)
      {
        case 0:
        case 5: value = 0; break;

        case 1:
        case 2: value = 1; break;

        case 3:
        case 4: value = 2; break;
      }

      break;

    case 2:
    case 3:

      switch (zm)
      {
        case 0:
        case 1: value = 2; break;

        case 2:
        case 3: value = 0; break;

        case 4:
        case 5: value = 1; break;
      }

      break;

    /* These cases are harder.  These blocks are divided diagonally
     * by the angled edges of the hexagons.  Some slope positive, and
     * others negative.  We flip the x value of the negatively sloped
     * pieces.  Then we check to see if the point in question falls
     * in the upper or lower half of the block.  That info, plus the
     * z status of the block determines the color.
     */

    case 1:
    case 4:

      /* Map the point into the block at the origin. */

      xl = xs-xm;
      zl = zs-zm;

      /* These blocks have negative slopes so we flip it horizontally. */

      if (((xm + zm) % 2) == 1)
      {
        xl = 1.0 - xl;
      }

      /* Avoid a divide-by-zero error. */

      if (xl == 0.0)
      {
        xl = 0.0001;
      }

      /* Is the angle less-than or greater-than 45 degrees? */

      brkindx = (zl / xl) < 1.0;

      /* was...
       * brkindx = (atan2(zl,xl) < (45 * M_PI_180));
       * ...but because of the mapping, it's easier and cheaper,
       * CPU-wise, to just use a good ol' slope.
       */

      switch (brkindx)
      {
        case TRUE:

          switch (zm)
          {
            case 0:
            case 3: value = 0; break;

            case 2:
            case 5: value = 1; break;

            case 1:
            case 4: value = 2; break;
          }

          break;

        case FALSE:

          switch (zm)
          {
            case 0:
            case 3: value = 2; break;

            case 2:
            case 5: value = 0; break;

            case 1:
            case 4: value = 1; break;
          }

          break;
      }
  }

  value = fmod(value, 3.0);

  return(value);
}



/*****************************************************************************
*
* FUNCTION
*
*   PickInCube(tv, p1)
*
* INPUT
*
*   ?
*
* OUTPUT
*   
* RETURNS
*
*   long integer hash function used, to speed up cacheing.
*   
* AUTHOR
*
*   Jim McElhiney
*   
* DESCRIPTION
*
*   A subroutine to go with crackle.
*
*   Pick a random point in the same unit-sized cube as tv, in a
*   predictable way, so that when called again with another point in
*   the same unit cube, p1 is picked to be the same.
*
* CHANGES
*
******************************************************************************/

static long PickInCube(tv, p1)
VECTOR tv, p1;
{
  int seed;
  VECTOR flo;

  /*
   * This uses floor() not FLOOR, so it will not be a mirror
   * image about zero in the range -1.0 to 1.0. The viewer
   * won't see an artefact around the origin.
   */

  flo[X] = floor(tv[X] - EPSILON);
  flo[Y] = floor(tv[Y] - EPSILON);
  flo[Z] = floor(tv[Z] - EPSILON);

  seed = Hash3d((int)flo[X], (int)flo[Y], (int)flo[Z]);

  POV_SRAND(seed);

  p1[X] = flo[X] + FRAND();
  p1[Y] = flo[Y] + FRAND();
  p1[Z] = flo[Z] + FRAND();

  return((long)seed);
}



/*****************************************************************************
*
* FUNCTION
*
*   Evaluate_Pattern
*
* INPUT
*
*   EPoint -- The point in 3d space at which the pattern
*   is evaluated.
*   TPat   -- Texture pattern struct
*   
* OUTPUT
*   
* RETURNS
*
*   DBL result usual 0.0 to 1.0 but may be 2.0 in hexagon
*   
* AUTHOR
*
*   adapted from Add_Pigment by Chris Young
*   
* DESCRIPTION
*
* CHANGES
*
******************************************************************************/

DBL Evaluate_TPat (TPat, EPoint)
TPATTERN *TPat;
VECTOR EPoint;
{
  DBL value = 0.0;
  VECTOR TPoint;

  Warp_EPoint (TPoint, EPoint, TPat);

  switch (TPat->Type)
  {
    case AGATE_PATTERN:    value = agate    (TPoint, TPat);   break;

    case BOZO_PATTERN:
    case SPOTTED_PATTERN:
    case BUMPS_PATTERN:    value = Noise    (TPoint);         break;

    case BRICK_PATTERN:    value = brick    (TPoint, TPat);   break;
    case CHECKER_PATTERN:  value = checker  (TPoint);         break;
    case CRACKLE_PATTERN:  value = crackle  (TPoint);         break;
    case GRADIENT_PATTERN: value = gradient (TPoint, TPat);   break;
    case GRANITE_PATTERN:  value = granite  (TPoint);         break;
    case HEXAGON_PATTERN:  value = hexagon  (TPoint);         break;
    case LEOPARD_PATTERN:  value = leopard  (TPoint);         break;
    case MANDEL_PATTERN:   value = mandel   (TPoint, TPat);   break;
    case MARBLE_PATTERN:   value = marble   (TPoint, TPat);   break;
    case ONION_PATTERN:    value = onion    (TPoint);         break;
    case RADIAL_PATTERN:   value = radial   (TPoint);         break;
    case SPIRAL1_PATTERN:  value = spiral1  (TPoint, TPat);   break;
    case SPIRAL2_PATTERN:  value = spiral2  (TPoint, TPat);   break;
    case PATTERN1_PATTERN: value = pattern1    (TPoint, TPat);   break;
    case PATTERN2_PATTERN: value = pattern2    (TPoint, TPat);   break;
    case PATTERN3_PATTERN: value = pattern3    (TPoint, TPat);   break;
    case WOOD_PATTERN:     value = wood     (TPoint, TPat);   break;

    case WAVES_PATTERN:    value = waves_pigm    (TPoint, TPat);   break;
    case RIPPLES_PATTERN:  value = ripples_pigm  (TPoint, TPat);   break;
    case WRINKLES_PATTERN: value = wrinkles_pigm (TPoint);   break;
    case DENTS_PATTERN:    value = dents_pigm    (TPoint);   break;
    case QUILTED_PATTERN:  value = quilted_pigm  (TPoint, TPat);   break;
    default: Error("Problem in Evaluate_TPat.");
  }

  if (TPat->Frequency !=0.0)
  {
    value = fmod(value * TPat->Frequency + TPat->Phase, 1.00001);
  }

  /* allow negative Frequency */

  if (value < 0.0)
  {
    value -= floor(value);
  }

  switch (TPat->Wave_Type)
  {
    case RAMP_WAVE:
      break;

    case SINE_WAVE:
      value = (1.0+cycloidal(value))*0.5;
      break;

    case TRIANGLE_WAVE:
      value = Triangle_Wave(value);
      break;

    case SCALLOP_WAVE:
      value = fabs(cycloidal(value*0.5));
      break;

    default: Error("Unknown Wave Type %d.",TPat->Wave_Type);
   }

  return(value);
}


/*****************************************************************************
*
* FUNCTION
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* DESCRIPTION
*
* CHANGES
*
******************************************************************************/

void Init_TPat_Fields (Tpat)
TPATTERN *Tpat;
{
  Tpat->Type       = NO_PATTERN;
  Tpat->Wave_Type  = RAMP_WAVE;
  Tpat->Flags      = NO_FLAGS;
  Tpat->References = 1;
  Tpat->Frequency  = 1.0;
  Tpat->Phase      = 0.0;
  Tpat->Warps      = NULL;
  Tpat->Next       = NULL;
  Tpat->Blend_Map  = NULL;
}



/*****************************************************************************
*
* FUNCTION
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* DESCRIPTION
*
* CHANGES
*
******************************************************************************/

void Copy_TPat_Fields (New, Old)
TPATTERN *New, *Old;
{
  *New = *Old;
  
  /* Copy warp chain */
  New->Warps = Copy_Warps(Old->Warps);

  New->Blend_Map = Copy_Blend_Map(Old->Blend_Map);

  /* Note, cannot copy Old->Next because we don't know what kind of
     thing this is.  It must be copied by Copy_Pigment, Copy_Tnormal etc.
  */

  if (Old->Type == BITMAP_PATTERN)
  {
     New->Vals.Image = Copy_Image(Old->Vals.Image);
  }
}



/*****************************************************************************
*
* FUNCTION
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* DESCRIPTION
*
* CHANGES
*
******************************************************************************/

void Destroy_TPat_Fields(Tpat)
TPATTERN *Tpat;
{
  Destroy_Warps(Tpat->Warps);
  Destroy_Blend_Map(Tpat->Blend_Map);
  /* Note, cannot destroy Tpat->Next nor pattern itself because we don't
     know what kind of thing this is.  It must be destroied by Destroy_Pigment, etc.
  */

  if (Tpat->Type == BITMAP_PATTERN)
  {
     Destroy_Image(Tpat->Vals.Image);
  }
}




/*****************************************************************************
*
* FUNCTION
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* DESCRIPTION
*
* CHANGES
*
******************************************************************************/

TURB *Create_Turb()
{
  TURB *New;

  New = POV_MALLOC(sizeof(TURB),"turbulence struct");

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

  New->Octaves = 6;
  New->Omega = 0.5;
  New->Lambda = 2.0;

  return(New);
}



/*****************************************************************************
*
* FUNCTION
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* DESCRIPTION
*
* CHANGES
*
******************************************************************************/

#if 0   /* Unused function [AED] */
static TURB *Copy_Turb(Old)
TURB *Old;
{
  TURB *New;

  if (Old != NULL)
  {
    New = Create_Turb();

    *New = *Old;
  }
  else
  {
    New=NULL;
  }

  return(New);
}
#endif


/*****************************************************************************
*
* FUNCTION
*
*   Translate_Tpattern
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   POV-Ray Team
*
* DESCRIPTION
*
* CHANGES
*
******************************************************************************/

void Translate_Tpattern(Tpattern,Vector)
TPATTERN *Tpattern;
VECTOR Vector;
{
  TRANSFORM Trans;

  if (Tpattern != NULL)
  {
    Compute_Translation_Transform (&Trans, Vector);

    Transform_Tpattern (Tpattern, &Trans);
  }
}



/*****************************************************************************
*
* FUNCTION
*
*   Rotate_Tpattern
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   POV-Ray Team
*
* DESCRIPTION
*
* CHANGES
*
******************************************************************************/

void Rotate_Tpattern(Tpattern,Vector)
TPATTERN *Tpattern;
VECTOR Vector;
{
  TRANSFORM Trans;

  if (Tpattern != NULL)
  {
    Compute_Rotation_Transform (&Trans, Vector);

    Transform_Tpattern (Tpattern, &Trans);
  }
}



/*****************************************************************************
*
* FUNCTION
*
*   Scale_Tpattern
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   POV-Ray Team
*
* DESCRIPTION
*
* CHANGES
*
******************************************************************************/

void Scale_Tpattern(Tpattern,Vector)
TPATTERN *Tpattern;
VECTOR Vector;
{
  TRANSFORM Trans;

  if (Tpattern != NULL)
  {
    Compute_Scaling_Transform (&Trans, Vector);

    Transform_Tpattern (Tpattern, &Trans);
  }
}



/*****************************************************************************
*
* FUNCTION
*
*   Transform_Tpattern
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   POV-Ray Team
*
* DESCRIPTION
*
* CHANGES
*
******************************************************************************/

void Transform_Tpattern(Tpattern,Trans)
TPATTERN *Tpattern;
TRANSFORM *Trans;
{
  WARP *Temp;

  if (Tpattern != NULL)
  {
    if (Tpattern->Warps == NULL)
    {
      Tpattern->Warps=Create_Warp(TRANSFORM_WARP);
    }
    else
    {
      if (Tpattern->Warps->Warp_Type != TRANSFORM_WARP)
      {
        Temp=Tpattern->Warps;

        Tpattern->Warps=Create_Warp(TRANSFORM_WARP);

        Tpattern->Warps->Next_Warp=Temp;
      }
    }

    Compose_Transforms (&( ((TRANS *)(Tpattern->Warps))->Trans), Trans);
  }
}


/*****************************************************************************
*
* FUNCTION
*
*   ripples_pigm
*
* INPUT
*
*   EPoint -- The point in 3d space at which the pattern
*   is evaluated.
*   TPat   -- Texture pattern struct
*
* OUTPUT
*   
* RETURNS
*
*   DBL value in the range 0.0 to 1.0
*   
* AUTHOR
*
*   POV-Ray Team
*
* DESCRIPTION   : Note this pattern is only used for pigments and textures.
*                 Normals have a specialized pattern for this.
*
* CHANGES
*   Nov 1994 : adapted from normal by [CY]
*
******************************************************************************/

static DBL ripples_pigm (EPoint, TPat)
VECTOR EPoint;
TPATTERN *TPat;
{
  register unsigned int i;
  register DBL length, index;
  DBL scalar =0.0;
  VECTOR point;

  for (i = 0 ; i < Number_Of_Waves ; i++)
  {
    VSub (point, EPoint, Wave_Sources[i]);
    VLength (length, point);

    if (length == 0.0)
      length = 1.0;

    index = length * TPat->Frequency + TPat->Phase;

    scalar += cycloidal(index);
  }

  scalar = 0.5*(1.0+(scalar / (DBL)Number_Of_Waves));

  return(scalar);
}


/*****************************************************************************
*
* FUNCTION
*
*   waves_pigm
*
* INPUT
*
*   EPoint -- The point in 3d space at which the pattern
*   is evaluated.
*   TPat   -- Texture pattern struct
*
* OUTPUT
*
* RETURNS
*
*   DBL value in the range 0.0 to 1.0
*
* AUTHOR
*
*   POV-Ray Team
*
* DESCRIPTION   : Note this pattern is only used for pigments and textures.
*                 Normals have a specialized pattern for this.
*
* CHANGES
*   Nov 1994 : adapted from normal by [CY]
*
******************************************************************************/

static DBL waves_pigm (EPoint, TPat)
VECTOR EPoint;
TPATTERN *TPat;
{
  register unsigned int i;
  register DBL length, index;
  DBL scalar = 0.0;
  VECTOR point;

  for (i = 0 ; i < Number_Of_Waves ; i++)
  {
    VSub (point, EPoint, Wave_Sources[i]);
    VLength (length, point);

    if (length == 0.0)
    {
      length = 1.0;
    }

    index = length * TPat->Frequency * frequency[i] + TPat->Phase;

    scalar += cycloidal(index)/frequency[i];
  }

  scalar = 0.2*(2.5+(scalar / (DBL)Number_Of_Waves));

  return(scalar);
}


/*****************************************************************************
*
* FUNCTION
*
*   dents_pigm
*
* INPUT
*
*   EPoint -- The point in 3d space at which the pattern
*   is evaluated.
*
* OUTPUT
*
* RETURNS
*
*   DBL value in the range 0.0 to 1.0
*   
* AUTHOR
*
*   POV-Ray Team
*   
* DESCRIPTION   : Note this pattern is only used for pigments and textures.
*                 Normals have a specialized pattern for this.
*
* CHANGES
*   Nov 1994 : adapted from normal by [CY]
*
******************************************************************************/

static DBL dents_pigm (EPoint)
VECTOR EPoint;
{
  DBL noise;

  noise = Noise (EPoint);

  return(noise * noise * noise);
}



/*****************************************************************************
*
* FUNCTION
*
*   wrinkles_pigm
*
* INPUT
*
*   EPoint -- The point in 3d space at which the pattern
*   is evaluated.
*   
* OUTPUT
*   
* RETURNS
*
*   DBL value in the range 0.0 to 1.0
*   
* AUTHOR
*
*   POV-Ray Team
*   
* DESCRIPTION   : Note this pattern is only used for pigments and textures.
*                 Normals have a specialized pattern for this.
*
* CHANGES
*   Nov 1994 : adapted from normal by [CY]
*
******************************************************************************/


static DBL wrinkles_pigm (EPoint)
VECTOR EPoint;
{
  register int i;
  DBL lambda = 2.0;
  DBL omega = 0.5;
  DBL value;
  VECTOR temp;

  value = Noise(EPoint);

  for (i = 1; i < 10; i++)
  {
    VScale(temp,EPoint,lambda);

    value += omega * Noise(temp);

    lambda *= 2.0;

    omega *= 0.5;
  }

  return(value/2.0);
}



/*****************************************************************************
*
* FUNCTION
*
*   quilted_pigm
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Dan Farmer & Chris Young
*   
* DESCRIPTION
*
* CHANGES
*
******************************************************************************/

static DBL quilted_pigm (EPoint, TPat)
VECTOR EPoint;
TPATTERN *TPat;
{
  VECTOR value;
  DBL t;

  value[X] = EPoint[X]-FLOOR(EPoint[X])-0.5;
  value[Y] = EPoint[Y]-FLOOR(EPoint[Y])-0.5;
  value[Z] = EPoint[Z]-FLOOR(EPoint[Z])-0.5;

  t = sqrt(value[X]*value[X]+value[Y]*value[Y]+value[Z]*value[Z]);

  t = quilt_cubic(t, TPat->Vals.Quilted.Control0, TPat->Vals.Quilted.Control1);

  value[X] *= t;
  value[Y] *= t;
  value[Z] *= t;

  return((fabs(value[X])+fabs(value[Y])+fabs(value[Z]))/3.0);
}



/*****************************************************************************
*
* FUNCTION
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* DESCRIPTION
*
* CHANGES
*
******************************************************************************/

#define INV_SQRT_3_4 1.154700538
DBL quilt_cubic(t,p1,p2)
DBL t,p1,p2;
{
 DBL it=(1-t);
 DBL itsqrd=it*it;
/* DBL itcubed=it*itsqrd; */
 DBL tsqrd=t*t;
 DBL tcubed=t*tsqrd;
 DBL val;

/* Originally coded as...

 val= (DBL)(itcubed*n1+(tcubed)*n2+3*t*(itsqrd)*p1+3*(tsqrd)*(it)*p2);

 re-written by CEY to optimise because n1=0 n2=1 always.
 
*/

 val = (tcubed + 3.0*t*itsqrd*p1 + 3.0*tsqrd*it*p2) * INV_SQRT_3_4;
 
 return(val);
}



/*****************************************************************************
*
* FUNCTION
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* DESCRIPTION
*
* CHANGES
*
******************************************************************************/

void Search_Blend_Map (value,Blend_Map,Prev,Cur)
DBL value;
BLEND_MAP *Blend_Map;
BLEND_MAP_ENTRY **Prev, **Cur;
{
  BLEND_MAP_ENTRY *P, *C;
  int Max_Ent=Blend_Map->Number_Of_Entries-1;

  /* if greater than last, use last. */

  if (value >= Blend_Map->Blend_Map_Entries[Max_Ent].value)
  {
    P = C = &(Blend_Map->Blend_Map_Entries[Max_Ent]);
  }
  else
  {
    P = C = &(Blend_Map->Blend_Map_Entries[0]);

    while (value > C->value)
    {
      P = C++;
    }
  }

  if (value == C->value)
  {
    P = C;
  }

  *Prev = P;
  *Cur  = C;
}



/*****************************************************************************
*
* FUNCTION
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* DESCRIPTION
*
* CHANGES
*
******************************************************************************/

static TURB *Search_For_Turb(Warps)
WARP *Warps;
{
  WARP* Temp=Warps;

  if (Temp!=NULL)
  {
    while (Temp->Next_Warp != NULL)
    {
      Temp=Temp->Next_Warp;
    }

    if (Temp->Warp_Type != CLASSIC_TURB_WARP)
    {
       Temp=NULL;
    }
  }

  return ((TURB *)Temp);
}