/*

  Ce programme genere des roches fractales
  pour POV Ray (en mesh)


  (c) 1998 Steven Pigeon

  modele de base: icosaedre (20 faces)

*/


#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

#include "utils.h"
#include "imath.h"



/******************************************************/
/*                                                    */
/* Contains a point                                   */
/*                                                    */
/******************************************************/
class Point
         {
          public:
           float alpha,beta,rayon;

           int operator==(Point &b);
           int operator>(Point &b);
           int operator<(Point &b);

           Point(); // constructeur par defaut
           Point(float a, float b, float r); // constructeur
           Point(Point & a); // copy-constructor
         };

/******************************************************/
Point Canonic(Point & a)
 {
  Point b;
  b.alpha=fmod(360.0f+a.alpha,360.0f);
  b.beta =fmod(360.0f+a.beta,360.0f);
  b.rayon=a.rayon;
  return b;
 }

/******************************************************/
Point::Point()
 {
  alpha=-1.0f;
  beta = 0.0f;
  rayon= 0.0f;
 }

/******************************************************/
Point::Point(float a, float b, float r)
 {
  alpha=a;
  beta =b;
  rayon=r;
 }

/******************************************************/
Point::Point(Point & a)
 {
  alpha=a.alpha;
  beta =a.beta;
  rayon=a.rayon;
 }

/******************************************************/
int Point::operator==(Point &b)
 {
  return (alpha==b.alpha) && (beta==b.beta);
 }

/******************************************************/
int Point::operator>(Point &b)
 {
  return (alpha>b.alpha) || ((alpha==b.alpha) && (beta>b.beta));
 }

/******************************************************/
int Point::operator<(Point &b)
 {
  return (alpha<b.alpha) || ((alpha==b.alpha) && (beta<b.beta));
 }

/******************************************************/
/*                                                    */
/* Global variables.                                  */
/*                                                    */
/******************************************************/
int MaxProf;
int mode=0;   // sharp or smoothed triangles ?

/******************************************************/
/*                                                    */
/* Look-up table functions and variables.             */
/*                                                    */
/******************************************************/
unsigned int   MaxPoint=0;        // sizeof table
unsigned int   PointPtr=0;       // pointer to next entry past last
Point         *PointTable=NULL; // point table

/******************************************************/
unsigned int PrimaryHash(Point & a)
 {
  unsigned int p1,p2;
  p1=*((unsigned int *)&a.alpha);
  p2=*((unsigned int *)&a.beta);
  return (p1 ^ p2 ) % MaxPoint; 
 }

/******************************************************/
Point Empty(-1.0f,0.0f,0.0f);

/******************************************************/
unsigned int PosTable(Point & a)
 {
  unsigned int i=PrimaryHash(a);

  while ( (!(PointTable[i]==Empty)) &&
          (!(PointTable[i]==a)))
   i = (i+21341u) % MaxPoint;

  return i;
 }

/******************************************************/
/*                                                    */
/* Checks if a point is in the look-up table          */
/*                                                    */
/******************************************************/
unsigned int PointExiste(Point & a)
 {
  if (PointPtr)
       {
        Point b=Canonic(a);
        unsigned int t=PosTable(b);
        return (PointTable[t]==b);
       }
  else return 0;
 }

/******************************************************/
/*                                                    */
/* Puts a point in the look-up table                  */
/*                                                    */
/******************************************************/
void PutPoint(Point & a)
 {
  Point b=Canonic(a);
  unsigned int t=PosTable(b);
  PointTable[t]=b;
  PointPtr++;
 }

/******************************************************/
/*                                                    */
/* Gets a point from the look-up table                */
/*                                                    */
/******************************************************/
void GetPoint(Point & a)
 {
  Point b=Canonic(a);
  unsigned int t=PosTable(b);
  a.rayon=PointTable[t].rayon;
 }

/******************************************************/
/*                                                    */
/* From spherical to rectangular coordinates          */
/*                                                    */
/******************************************************/
float Polar2Cart_X(Point & a)
 {
  return (float)(a.rayon * sin(a.beta/57.2957) * cos(a.alpha/57.2957));
 }

/******************************************************/
float Polar2Cart_Y(Point & a)
 {
  return (float)(a.rayon * cos(a.beta/57.2957));
 }

/******************************************************/
float Polar2Cart_Z(Point & a)
 {
  return (float)(a.rayon * sin(a.alpha/57.2957) * sin(a.beta/57.2957));
 }

/******************************************************/
Point Normalise(float x, float y, float z)
 {
  Point a;
  float norm = sqrt(x*x+y*y+z*z);

  a.alpha = x/norm;
  a.beta  = y/norm;
  a.rayon = z/norm;

  return a;
 }

/******************************************************/
void InstancieTriangle(Point & a, Point & b, Point & c)
 {
  if (mode==0)
       printf("triangle{<%.5f,%.5f,%.5f>,<%.5f,%.5f,%.5f>,<%.5f,%.5f,%.5f>}\n",
              Polar2Cart_X(a),Polar2Cart_Y(a),Polar2Cart_Z(a),
              Polar2Cart_X(b),Polar2Cart_Y(b),Polar2Cart_Z(b),
              Polar2Cart_X(c),Polar2Cart_Y(c),Polar2Cart_Z(c) );
  else {
        Point N1 = Normalise(Polar2Cart_X(a),Polar2Cart_Y(a),Polar2Cart_Z(a));
        Point N2 = Normalise(Polar2Cart_X(b),Polar2Cart_Y(b),Polar2Cart_Z(b));
        Point N3 = Normalise(Polar2Cart_X(c),Polar2Cart_Y(c),Polar2Cart_Z(c));
        printf("smooth_triangle{<%.4f,%.4f,%.4f>,<%.2f,%.2f,%.2f>,"
               "<%.4f,%.4f,%.4f>,<%.2f,%.2f,%.2f>,"
               "<%.4f,%.4f,%.4f>,<%.2f,%.2f,%.2f>}\n",
                Polar2Cart_X(a),Polar2Cart_Y(a),Polar2Cart_Z(a),N1.alpha,N1.beta,N1.rayon,
                Polar2Cart_X(b),Polar2Cart_Y(b),Polar2Cart_Z(b),N2.alpha,N2.beta,N2.rayon,
                Polar2Cart_X(c),Polar2Cart_Y(c),Polar2Cart_Z(c),N3.alpha,N3.beta,N3.rayon
               );
       }
 }

float Expos=1.0f;
/******************************************************/
float Perturbation(int profondeur)
 {
  return ( (rand()/(float)RAND_MAX) - (rand()/(float)RAND_MAX) ) /
           exp( Expos*(1+MaxProf-profondeur) );
 }

/******************************************************/
float Perturbe(Point & a, Point & b, int profondeur)
 {
  return ((a.rayon+b.rayon)/2.0f) + Perturbation(profondeur);
 }

/******************************************************/
void Recurse(Point & a, Point & b, Point & c, int profondeur, int Polar=0)
 {
  if (profondeur)
      { 
       // recursons joyeusement
       Point d,e,f;

       if (Polar)
         // north polar region
        {
         d.alpha=b.alpha;
         d.beta=(a.beta+b.beta)/2;
         if (PointExiste(d))
              GetPoint(d);
         else {
               d.rayon=Perturbe(a,b,profondeur);
               PutPoint(d);
              }

         e.alpha=c.alpha;
         e.beta=(a.beta+c.beta)/2;
         if (PointExiste(e)) 
              GetPoint(e);
         else {
               e.rayon=Perturbe(a,c,profondeur);
               PutPoint(e);
              }
         
         f=Point( (b.alpha+c.alpha)/2.0f, (b.beta+c.beta)/2.0f, 0.0f );
         if (PointExiste(f))
              GetPoint(f);
         else {
               f.rayon=Perturbe(b,c,profondeur);
               PutPoint(f);
              }

         Recurse(a,d,e,profondeur-1,1);
         Recurse(d,e,f,profondeur-1);
         Recurse(d,b,f,profondeur-1);
         Recurse(e,f,c,profondeur-1);
        }
       else
         {
          d=Point( (a.alpha+b.alpha)/2.0f, (a.beta+b.beta)/2.0f, 0.0f );
          if ( PointExiste(d))
               GetPoint(d);
          else {
                d.rayon=Perturbe(a,b,profondeur);
                PutPoint(d);
               }

          e=Point( (a.alpha+c.alpha)/2.0f,(a.beta+c.beta)/2.0f,0.0f);
          if ( PointExiste(e))
               GetPoint(e);
          else {
                e.rayon=Perturbe(a,c,profondeur);
                PutPoint(e);
               }

          f=Point( (b.alpha+c.alpha)/2.0f,(b.beta+c.beta)/2.0f,0.0f);
          if ( PointExiste(f))
               GetPoint(f);
          else {
                f.rayon=Perturbe(b,c,profondeur);
                PutPoint(f);
               }

          Recurse(a,d,e,profondeur-1);
          Recurse(d,e,f,profondeur-1);
          Recurse(d,b,f,profondeur-1);
          Recurse(e,f,c,profondeur-1);
         }
      }
  else InstancieTriangle(a,b,c);
 }

/******************************************************/
void main( int argc, char *argv[])
 {

  if (argc<2)
        printf(" Usage:\n\n"
               " PovRockGen <depth> <randseed> <smoothness> <mode>\n\n"
               " <depth>      : recursion depth. Legal values: 0..8\n"
               " <randseed>   : either a number (the seed) or '+' to randomize\n"
               " <smoothness> : a real number, negative: 0.0 is more choatic, \n"
               "                the biggest, the more smooth. Defaults to 1.0. \n"
               " <mode>       : if '+', uses smoothed triangles.\n\n"
               " Redirect stdout to a file. For e.g.,\n"
               " PovRockGen 6 + 0.5 + > z.pov\n\n");
  else {
        if (argc>1) 
             MaxProf=max(0,min(8,atoi(argv[1])));
        else MaxProf=3;

        if (argc>2)
           {
            if  (argv[2][0]=='+') 
                 {
                  int i = time(NULL);
                  fprintf(stderr,"seed:%d\n",i);
                  srand(i);
                 }
            else srand(atoi(argv[2]));
           }
      
        if (argc>3) Expos=max(0.0f,(float)atof(argv[3]));

        if (argc>4) mode = (argv[4][0]=='+');
  

        fprintf(stderr,"mesh has %d triangles\n", 20 * ExpoDiscrete(4,MaxProf));
        MaxPoint = 10 * ExpoDiscrete(2,MaxProf)*(ExpoDiscrete(2,MaxProf)+1);
        fprintf(stderr,"mesh has ~%d summits\n",MaxPoint);

        PointTable = new  Point[MaxPoint];
        PointPtr=0;

        printf("mesh\n {\n");

          // Recurse from an initial shape: an icosahedron

          // couronne
          Recurse(Point(0,0,1),Point(-36,60,1),Point( +36,60,1), MaxProf,1);  fprintf(stderr,".");
          Recurse(Point(72,0,1),Point(+36,60,1),Point( 108,60,1), MaxProf,1); fprintf(stderr,".");
          Recurse(Point(144,0,1),Point(108,60,1),Point(180,60,1), MaxProf,1); fprintf(stderr,".");
          Recurse(Point(-72,0,1),Point(-108,60,1),Point(-36,60,1), MaxProf,1); fprintf(stderr,".");
          Recurse(Point(-144,0,1),Point(-180,60,1),Point(-108,60,1), MaxProf,1); fprintf(stderr,".");
  
          // ceinture
          Recurse(Point(0,120,1),Point(-36,60,1),Point( +36,60,1), MaxProf);  fprintf(stderr,".");
          Recurse(Point(72,120,1),Point(+36,60,1),Point( 108,60,1), MaxProf); fprintf(stderr,".");
          Recurse(Point(144,120,1),Point(108,60,1),Point(180,60,1), MaxProf); fprintf(stderr,".");
          Recurse(Point(-72,120,1),Point(-108,60,1),Point(-36,60,1), MaxProf); fprintf(stderr,".");
          Recurse(Point(-144,120,1),Point(-180,60,1),Point(-108,60,1), MaxProf); fprintf(stderr,".");

          Recurse(Point(36,60,1),Point(0,120,1),Point(72,120,1), MaxProf); fprintf(stderr,".");
          Recurse(Point(108,60,1),Point(72,120,1),Point(144,120,1), MaxProf); fprintf(stderr,".");
          Recurse(Point(-180,60,1),Point(-144,120,1),Point(-216,120,1), MaxProf); fprintf(stderr,".");
          Recurse(Point(-36,60,1),Point(0,120,1),Point(-72,120,1), MaxProf); fprintf(stderr,".");
          Recurse(Point(-108,60,1),Point(-72,120,1),Point(-144,120,1), MaxProf); fprintf(stderr,".");
  
          // base
          Recurse(Point(36,180,1),Point(0,120,1),Point(72,120,1), MaxProf,1); fprintf(stderr,".");
          Recurse(Point(108,180,1),Point(72,120,1),Point(144,120,1), MaxProf,1); fprintf(stderr,".");
          Recurse(Point(-180,180,1),Point(-144,120,1),Point(-216,120,1), MaxProf,1); fprintf(stderr,".");
          Recurse(Point(-36,180,1),Point(0,120,1),Point(-72,120,1), MaxProf,1); fprintf(stderr,".");
          Recurse(Point(-108,180,1),Point(-72,120,1),Point(-144,120,1), MaxProf,1); fprintf(stderr,".\n");

        printf(" }\n");

        delete [] PointTable;
       }
 }
