//////////////////////////////////////////////////////////////////////////////
//     File: hypercloud.h
//  Created: 28.11.2001 in Delphi for I-Tex 2.2, converted to C++ 3.2.2002
//   Author: Tobias Reichert
//    Email: tobias_reichert@i-tex.de
//  Summary: The class THyperCloudTexture is initialized
//           with the Init(...) method that sets random seed and
//           initializes the parameters and random value cache.
//           double GetValue4D(x,y,z,time) returns a depth value [0..1).
//           Animated this looks life boiling water ;-)
//           Computation Expense : O(Octaves*Density*2^Dimension),
//           nothing for realtime apps, but can be optimized a lot
//////////////////////////////////////////////////////////////////////////////

#include "math.hpp"
#include "math.h"

class THyperCloudTexture
      {
       int RSeed;           //Random Seed, use full range
       int Octaves;         //Detail, Proportional to Computation Expense, use 6+
       int Density;         //Number of Puffs per Gridcell, use [1..20]
       double SubScale;     //Resolution for Octave(n)=Subscale^n, use [1 .. 4]
       double Persistance;  //Valence for Octave(n)=
                            //Persistance^n/Sum(all Octave Valences)
                            //use [0 .. 1)
       double ReciPuffSize; //1/Size of Cloud Puffs, use [0.5 .. 1)
       double Amplify;      //Values can be very low, scale with this
       double Bias;         //Make also Brighter with this Value

       double ReciOctaveSum; //is 1/Sum(all Octave Valances), Initialized by Init()

       double FNoiseArray[1024]; //RandomPool of Doubles
       int INoiseArray[1024]; //RandomPool of Integers

       private:
        double __fastcall DirectRandom(int i);

       public:
        void __fastcall Init(int    pRSeed,
                             int    pOctaves,
                             int    pDensity,
                             double pSubScale,
                             double pPersistance,
                             double pPuffSize,
                             double pAmplify,
                             double pBias);

        double __fastcall GetValue4D(double x,
                                     double y,
                                     double z,
                                     double t);
      };


double __fastcall THyperCloudTexture::DirectRandom(int i)
{
 return abs(((i +  i * ( 89603  +   i *46703) ) * 114269) % 9995057 )
            * (double)1.000494544453323e-7;
}




void __fastcall THyperCloudTexture::Init(int    pRSeed,
                                         int    pOctaves,
                                         int    pDensity,
                                         double pSubScale,
                                         double pPersistance,
                                         double pPuffSize,
                                         double pAmplify,
                                         double pBias)
 //Example :  THyperCloudTexture *Clouds;
 //           Clouds=new(THyperCloudTexture);
 //           Clouds->Init(123,6,20,2,0.6,0.8,4,0.1);

{
 //Copy Parameters to Internal Variables

 RSeed = pRSeed;
 Octaves = pOctaves;
 Density = pDensity;
 SubScale = pSubScale;
 Persistance = pPersistance;

 //Puffs with Size=0 would not be very useful (and cause DbZ)
 if (pPuffSize>0) { ReciPuffSize = sqrt(1/pPuffSize); }
                  else
                  { ReciPuffSize = 10; }

 Amplify = pAmplify;
 Bias = pBias;

 //Compute 1/Sum(persistance^n,n=0..octaves-1);

 if (Octaves==0) {Octaves=1;} //Just for the case ...

 if (Persistance==1)
          { ReciOctaveSum = 1/Octaves; }
          else
          { ReciOctaveSum = (Persistance-1) / (IntPower(Persistance,Octaves) - 1) ; }

 //Init Cached Random Arrays

 for (int i=0; i<1024; i++)
 {
  FNoiseArray[i] = DirectRandom(i*101+RSeed);
  INoiseArray[i] = (int)(33111*DirectRandom(i*91+RSeed));
 }
}

double __fastcall THyperCloudTexture::GetValue4D(double x,
                                                 double y,
                                                 double z,
                                                 double t)
{
 int i,d; //Some Working Variables
 int tx,ty,tz,tt; //Scan loop Variables
 int nd; //Octave "ANDing" Argument = 2^n-1

 double c=0; //Return Value
 double b;
 double p=1; //Persistance

 int    intx,inty,intz,intt;
 int    offx,offy,offz,offt;
 double fracx,fracy,fracz,fract;
 double fx,fy,fz,ft;
 int inoff;
 double SqrDist;
 double PuffDens;
 double ReciSize;

 int cachex,cachey,cachez,cachet;
 int r;

 for (i=0;i<Octaves;i++) //Add Octaves
 {
  nd=(256<<i)-1;


  inoff = i * 44111 + RSeed; //Every Octave uses modified RSeed

  intx  = (int)x; //round down to Int of x,y,z,t
  inty  = (int)y;
  intz  = (int)z;
  intt  = (int)t;

  fracx = x-(double)intx; //fraction [0..1) value of x,y,z,t
  fracy = y-(double)inty;
  fracz = z-(double)intz;
  fract = t-(double)intt;

  //Determine Cells to be checked
  (fracx<0.5) ? (offx=-1) : (offx=0);
  (fracy<0.5) ? (offy=-1) : (offy=0);
  (fracz<0.5) ? (offz=-1) : (offz=0);
  (fract<0.5) ? (offt=-1) : (offt=0);

  //Initialize Octave-Value b
  b=0;

  for (d=0;d<Density;d++)
  {

   for (tt=offt;tt<=offt+1;tt++) //Scan for Overlapping Cells in dim-(t)ime
   {
    cachet = ((intt+tt) & nd);

    for (tz=offz;tz<=offz+1;tz++) //Scan for Overlapping Cells in dim-z
    {
     cachez = ((intz+tz) & nd);

     for (ty=offy;ty<=offy+1;ty++) //Scan for Overlapping Cells in dim-y
     {
      cachey = ((inty+ty) & nd);

      for (tx=offx;tx<=offx+1;tx++) //Scan for Overlapping Cells in dim-x
      {
       cachex = ((intx+tx) & nd);

       //Compute x-coordinate of element d
       r = cachex*33001 + cachey*49937 + cachez*53441 + cachet*21839 + inoff;
       fx = (double)tx+FNoiseArray[(INoiseArray[r & 1023] + r * 52223) & 1023]-fracx;

       //Compute y-coordinate of element d
       r = cachex*49937 + cachey*53441 + cachez*21839 + cachet*33001 + inoff;
       fy = (double)ty+FNoiseArray[(INoiseArray[r & 1023] + r * 21839) & 1023]-fracy;

       //Compute t-coordinate of element d
       r = cachex*53441 + cachey*21839 + cachez*33001 + cachet*49937 + inoff;
       fz = (double)tz+FNoiseArray[(INoiseArray[r & 1023] + r * 49937) & 1023]-fracz;

       //Compute t-coordinate of element d
       r = cachex*21839 + cachey*33001 + cachez*49937 + cachet*53441 + inoff;
       ft = (double)tt+FNoiseArray[(INoiseArray[r & 1023] + r * 53441) & 1023]-fract;

       //Compute 1/Size of element d
       r = cachex*33001 + cachey*49937 + cachez*53441 + cachet*21839 + inoff;
       ReciSize = ((2+FNoiseArray[(INoiseArray[r & 1023] + r * 52223) & 1023])*ReciPuffSize);

       //Compute brightness of element d
       r = cachex*33001 + cachey*21839 + cachez*53441 + cachet*49937 + inoff;
       PuffDens = 0.1+0.9*(FNoiseArray[(INoiseArray[r & 1023] + r * 52223) & 1023]);

       //Compute 1-(Pixel-Puff-Distance)^2 , Distance = Sum((Dimension-Distance)^2)
       SqrDist = PuffDens * (1-(fx*fx + fy*fy + fz*fz + ft*ft)*ReciSize*ReciSize);
       //If distance is lower (1-distance greater) than beforce, remember it
       if (SqrDist>b) { b=SqrDist; }

      } //End tx-loop
     } //End ty-loop
    } //End tz-loop
   } //End tt-loop

   //Change Random for each Element
   inoff += 101;

  } //End d-loop

  //Add Result of the Octave * Persistance
  c += sqrt(b)*p;

  //Higher Octaves are scaled with the SubScale Parameter
  x*=SubScale;
  y*=SubScale;
  z*=SubScale;
  t*=SubScale;

  //Compute Next Persistance
  p*=Persistance;
 }

 //Scale and Transform Result and Return Value
 return  1-exp(-c*Amplify * ReciOctaveSum - Bias);
}