/*
      Rendu de bidules freedirs,
               Par ChewBaccA !

   (c) personne

   message a kwa : si tu me lis, je suis desole pour ce code
                   immonde et sale :)
*/

#include <math.h>
#include <allegro.h>

#define FOV 120

// definition de types
typedef struct
{
  float x, y, z;
} VECTEUR;

typedef struct
{
  unsigned short x, y;
} point;

typedef float TMatrice [3][3];       // matrice 3x3

/////////////////////////////////////////////////////////////////////
//
// variables globales
//
/////////////////////////////////////////////////////////////////////

BITMAP *virscr;         // buffer dans lequel on ecrit (VIRtual SCReen)
PALETTE pal;

////////////////////////////////////////////////////////////////////
//
// position & direction de la camera
//
////////////////////////////////////////////////////////////////////

VECTEUR camera    = { 0,   0,   0 }; // position de la camera

TMatrice cam_matrix;                 // matrice de rotation de la camera

float cam_rotx = 0,
      cam_roty = 0,
      cam_rotz = 0;                  // angle de la rotation de la direction

////////////////////////////////////////////////////////////////////

unsigned char tex[65536];            // texture
point interpol[321][201];            // coordonnees a interpoler

// declarations de fonctions
void calc_camera_matrix();
void tunnel_free_dir();              // calcule un tunnel (cylindre)
void sphere_free_dir();              // calcule une sphere
void plane_free_dir();               // calcule deux plans

void interpolation();                // interpole le tableau interpol

void prod_vectmat(VECTEUR *v, VECTEUR *v1, TMatrice m1);
void rot_mat(TMatrice m, float rx, float ry, float rz);
void normalize(VECTEUR *v);          // ramene v a une longueur de 1

/////////////////////////////////////////////////////////////////////////////

int main() {
  int x, y;
  unsigned char c;

  float alpha=0;
  unsigned char speed=15;

  int freedir_type=1;            // type d'objet rendu
  int timbre_poste=1;

  // initialisation (on s'en branle de ca :) )
  allegro_init();
  install_keyboard();
  install_mouse();
  set_gfx_mode(GFX_VGA, 320, 200, 0, 0);

  virscr = create_bitmap(320, 200);
  clear(virscr);

  for(y=0;y<256;y++)
    for(x=0;x<256;x++)
      tex[(y<<8)+x] = x^y;    // on g‚nŠre une texture … 2 francs

  for(x=0;x<256;x++)
    pal[x].r = pal[x].g = pal[x].b = x>>2; // palette en niveaux de gris
  set_palette(pal);

  // c'est a partir de la que c'est interessant

  do {
    // modifier cam et direction
    camera.z += speed;
    camera.x = (80*cos(alpha/2)  + 60*sin(alpha/3));
    camera.y = (100*sin(alpha/4) + 40*cos(alpha));
    alpha += .26;

    cam_rotx += .013;    //   On change les angles des rotations
    cam_roty += .004;    //
    cam_rotz += .009;    //

    rot_mat(cam_matrix, cam_rotx, cam_roty, cam_rotz);
    //  et on recalcule la matrice de rotation de la cam

    // calculer l'espace tous les 8 pixels...
    switch(freedir_type)
    {
       case 1 : tunnel_free_dir();
                break;
       case 2 : sphere_free_dir();
                break;
       case 3 : plane_free_dir();
                break;
    }

    // ... puis on interpole ce qu'on a deja calcule
    interpolation();

    // on affiche le resultat
    for(y=0;y<200;y++)
      for(x=0;x<320;x++)
      {
        c = tex[((interpol[x][y].y<<8)+interpol[x][y].x)&0xffff];
        virscr->line[y][x] = c;
      }

    if(timbre_poste)
      for(y=0;y<200;y+=8)
        for(x=0;x<320;x+=8)
          _putpixel(virscr, x>>3, y>>3, tex[((interpol[x][y].y<<8)+interpol[x][y].x)&0xffff] );

    blit(virscr, screen, 0, 0, 0, 0, 320, 200);

    // gestion du clavier
    if(key[KEY_F1])
      freedir_type = 1;
    if(key[KEY_F2])
    {
      freedir_type = 2;
      speed=0;
      camera.x = camera.y = camera.z = 0;
    }
    if(key[KEY_F3])
      freedir_type = 3;

    if(key[KEY_BACKSPACE])
      timbre_poste = 1-timbre_poste;

    if(key[KEY_PLUS_PAD])
      speed++;
    if(key[KEY_MINUS_PAD])
      speed--;

  } while(!key[KEY_ESC]);

  allegro_exit();
  return 0;
}


void tunnel_free_dir()
// remplit interpol pour faire un tunnel
{
   int x, y;
   VECTEUR d;           // vecteur direction pour chaque rayon
   VECTEUR inter;       // point d'intersection entre le rayon et le tunnel

   float a, b, c;       // coefficient du polynome
   float delta;         // discriminant du polynome
   float t, t1, t2;     // t solutions

   static int rayon=256;// rayon du tunnel

   // calculs constants pour l'image
   c = (camera.x * camera.x) + (camera.y * camera.y) - (rayon*rayon);
   // c = (cxý+cyý-rý)

   for(y=0;y<201;y+=8)
     for(x=0;x<321;x+=8)
     // tous les 8 pixels
     {
       // on calcule la direction
       d.x = ((float)(x-160))/FOV;
       d.y = ((float)(y-100))/FOV;
       d.z = 1;
       prod_vectmat(&d, &d, cam_matrix);

       // on ramene la norme de d a 1
       normalize(&d);

       a = (d.x*d.x) + (d.y*d.y);
       // a = dxý+dyý
       b = 2*(camera.x*d.x + camera.y*d.y);
       // b = 2*(cx*dx+cy*dy)

       delta = (b*b) - 4*a*c;
       // classique : delta = bý-4ac
       if(delta>0) {
         delta = sqrt(delta);
         // on calcule racine(delta) avant pour economiser un calcul de racine

         t1 = (-b-delta)/(2*a);
         //       -b - racine(delta)
         // t1 = --------------------
         //            2a
         if(t1>0) t=t1;
         // t est la solution positive

         else
         {
           t2 = (-b+delta)/(2*a);
           //       -b + racine(delta)
           // t2 = --------------------
           //            2a
           t = t2;
         }

         // maintenant qu'on a t, on peut calculer les coordonnees du point
         // d'intersection entre le tunnel et le rayon
         // inter = camera + t*direction
         inter.x = camera.x + t*d.x;
         inter.y = camera.y + t*d.y;
         inter.z = camera.z + t*d.z;

         interpol[x][y].x =
           (unsigned short) (fabs( (256 * atan2(inter.y, inter.x) ) /PI) );
         interpol[x][y].y =
           (unsigned short)  fabs( inter.z );

       }
       else
       {
       // si il n'y a pas d'intersection
         interpol[x][y].x = 0;
         interpol[x][y].y = 0;
       }
     }
}

void plane_free_dir()
// remplit interpol pour faire deux plans
{
   int x, y;
   VECTEUR d;           // vecteur direction pour chaque rayon
   VECTEUR inter;       // point d'intersection entre le rayon et le tunnel

   float t;             // t solution

   static int hauteur=256;  // distance des plans sur Oy

   for(y=0;y<201;y+=8)
     for(x=0;x<321;x+=8)
     // tous les 8 pixels
     {
       // on calcule la direction
       d.x = ((float)(x-160))/FOV;
       d.y = ((float)(y-100))/FOV;
       d.z = 1;
       prod_vectmat(&d, &d, cam_matrix);

       // on ramene la norme de d a 1
       normalize(&d);

       if(d.y>0)
       // on regarde vers le plan du haut
       {
          t = (hauteur - camera.y) / d.y;

          // maintenant qu'on a t, on peut calculer les coordonnees du point
          // d'intersection entre le tunnel et le rayon
          // inter = camera + t*direction
          inter.x = camera.x + t*d.x;
          inter.y = camera.y + t*d.y;
          inter.z = camera.z + t*d.z;

          // mapping planaire
          interpol[x][y].x =
            (unsigned short) fabs( inter.x );
          interpol[x][y].y =
            (unsigned short) fabs( inter.z );
       }
       else if(d.y<0)
       // on regarde vers le plan du bas
       {
          t = (-hauteur - camera.y) / d.y;

          inter.x = camera.x + t*d.x;
          inter.y = camera.y + t*d.y;
          inter.z = camera.z + t*d.z;

          interpol[x][y].x =
            (unsigned short) fabs( inter.x );
          interpol[x][y].y =
            (unsigned short) fabs( inter.z );
       }
       else
       {
       // si il n'y a pas d'intersection
         interpol[x][y].x = 0;
         interpol[x][y].y = 0;
       }
     }
}

void sphere_free_dir()
// remplit interpol pour faire une sphere
{
   int x, y;
   VECTEUR d;           // vecteur direction pour chaque rayon
   VECTEUR inter;       // point d'intersection entre le rayon et le tunnel

   float a, b, c;       // coefficient du polynome
   float delta;         // discriminant du polynome
   float t, t1, t2;     // t solutions

   static int rayon=600;   // rayon de la sphere

   c = (camera.x * camera.x) + (camera.y * camera.y) +(camera.z * camera.z) - (rayon*rayon);
   // c = (cxý+cyý+czý-rý)

   for(y=0;y<201;y+=8)
     for(x=0;x<321;x+=8)
     // tous les 8 pixels
     {
       // on calcule la direction
       d.x = ((float)(x-160))/FOV;
       d.y = ((float)(y-100))/FOV;
       d.z = 1;
       prod_vectmat(&d, &d, cam_matrix);

       // on ramene la norme de d a 1
       normalize(&d);

       a = (d.x*d.x) + (d.y*d.y) + (d.z*d.z);
       // a = dxý+dyý+dzý
       b = 2*(camera.x*d.x + camera.y*d.y + camera.z*d.z);
       // b = 2*(cx*dx+cy*dy+cz*dz)

       delta = (b*b) - 4*a*c;
       // classique : delta = bý-4ac
       if(delta>0)
       {
         delta = sqrt(delta);
         // on calcule racine(delta) avant pour economiser un calcul de racine

         t1 = (-b-delta)/(2*a);
         //       -b - racine(delta)
         // t1 = --------------------
         //            2a

         t2 = (-b+delta)/(2*a);
         //       -b + racine(delta)
         // t2 = --------------------
         //            2a

         t = (t1>0)?t1 : t2;
         // t est la solution positive

         // maintenant qu'on a t, on peut calculer les coordonnees du point
         // d'intersection entre le tunnel et le rayon
         // inter = camera + t*direction
         inter.x = camera.x + t*d.x;
         inter.y = camera.y + t*d.y;
         inter.z = camera.z + t*d.z;

         // on determine les coordonnees (u,v) de la texture mappee sur le
         interpol[x][y].x =
           asin( (inter.y) / rayon )*128+128;
         interpol[x][y].y =
           (unsigned short) (fabs( (256 * atan2(inter.z, inter.x) ) /PI) );

       }
       else
       {
       // si il n'y a pas d'intersection
         interpol[x][y].x = 0;
         interpol[x][y].y = 0;
       }
     }
}

/////////////////////////////////////////////////////////////////////////
//
//  OPERATIONS SUR LES VECTEURS ...
//
/////////////////////////////////////////////////////////////////////////

void normalize(VECTEUR *v)
// ramene v a une longueur de 1
{
  float l;

  l = sqrt((v->x*v->x) + (v->y*v->y) + (v->z*v->z));
  v->x = v->x/(l?l:1);
  v->y = v->y/(l?l:1);
  v->z = v->z/(l?l:1);
}

void prod_vectmat(VECTEUR *v, VECTEUR *v1, TMatrice m1)
// multiplie v1 par m1 et place le resultat dans v
{
   int i;
   float temp[3] = { 0, 0, 0 };

   for(i=0;i<3;i++)
   {
     temp[i] += v1->x * m1[i][0];
     temp[i] += v1->y * m1[i][1];
     temp[i] += v1->z * m1[i][2];
   }

   v->x = temp[0];
   v->y = temp[1];
   v->z = temp[2];
}

void rot_mat(TMatrice m, float rx, float ry, float rz)
{
  float cosx, sinx;
  float cosy, siny;
  float cosz, sinz;

  cosx = (float) cos(rx);
  sinx = (float) sin(rx);
  cosy = (float) cos(ry);
  siny = (float) sin(ry);
  cosz = (float) cos(rz);
  sinz = (float) sin(rz);

  m[0][0] = cosy*cosz;
  m[0][1] = cosy*sinz;
  m[0][2] = -siny;

  m[1][0] = sinx*siny*cosz - cosx*sinz;
  m[1][1] = sinx*siny*sinz + cosx*cosz;
  m[1][2] = sinx*cosy;

  m[2][0] = cosx*siny*cosz + sinx*sinz;
  m[2][1] = cosx*siny*sinz - sinx*cosz;
  m[2][2] = cosx*cosy;
}

/////////////////////////////////////////////////////////////////////////////

void interpolation()
// interpole les coordonnees x et y de offset_x et offset_y dans virscr
// merci a cyg pour sa doc  (TROP MERCI Cyg ! :) )
{
/*
   formule d'interpolation :
     v(x,y) = ( (8-x)*(8-y)*c1 + x*(8-y)*c2 + (8-x)*y*c3 + x*y*c4 ) / 64

   d'ou :
     dv(x,y)/dx = ( -(8-y)*v1 + (8-y)*v2 - y*v3 + y*v4 ) / 64
   et
     dv(x,y)/dy = ( -(8-x)*v1 - x*v2 + (8-x)*v3 + x*v4 ) / 64

   dv/dx pas constant : il varie avec y, on calcule
                        donc la variation de dýv(x,y)/dxy

     dýv(x,y)/dxy = ( v1 - v2 - v3 + v4 ) / 64
*/

  int _x, _y;             // coordonnees dans le tableau de vals a interpoler
  int x, y;               // coordonnees a l'interieur du carre
  unsigned long  v;       // valeur courante (interpolee)
  unsigned long _v;       // valeur du debut de la ligne
  unsigned long dvdx;     // variation selon x (= dv(x,y)/dx )
  unsigned long dvdy;     // variation selon y (= dv(x,y)/dy )
  unsigned long dvdxdy;   // variation de dvdx quand y varie
  unsigned long v1, v2, v3, v4;   // valeurs aux 4 coins du carres

  // interpolation des x
  for(_y=0;_y<200;_y+=8)      // on scanne les differentes valeurs
    for(_x=0;_x<320;_x+=8)    // a interpoler
    // et on interpole les carres 8x8
    {
      v1 = interpol[_x][_y].x<<8;     // on recupere les valeurs aux 4 coins
      v2 = interpol[_x+8][_y].x<<8;   // du carre a interpoler
      v3 = interpol[_x][_y+8].x<<8;   // (on les passe sur 16 bits pour avoir
      v4 = interpol[_x+8][_y+8].x<<8; // plus de precision dans la variation)

      _v = v = v1;                 // on initialise le x a interpoler

      dvdx   = (-v1+v2)>>3;        // variation en x quand y=0
      dvdy   = (-v1+v3)>>3;        // variation en y quand x=0
      dvdxdy = (v1-v2-v3+v4)>>6;   // variation de dvdx quand y varie

      for(y=_y;y<_y+8;y++) {       // on parcourt le carre de haut en bas
        for(x=_x;x<_x+8;x++) {     // de gauche a droite
          v += dvdx;               // on fait varier la valeur a interpoler
          interpol[x][y].x = v>>8; // on repasse la valeur sur 8 bits
        }
        _v += dvdy;                // on fait varier la valeur de gauche
        v = _v;                    // on recommence au debut de la ligne
        dvdx += dvdxdy;            // on fait varier dvdx
      }

    }

  // interpolation des y
  for(_y=0;_y<200;_y+=8)      // on scanne les differentes valeurs
    for(_x=0;_x<320;_x+=8)    // a interpoler
    // et on interpole les carres 8x8
    {
      v1 = interpol[_x][_y].y<<8;     // on recupere les valeurs aux 4 coins
      v2 = interpol[_x+8][_y].y<<8;   // du carre a interpoler
      v3 = interpol[_x][_y+8].y<<8;   // (on les passe sur 16 bits pour avoir
      v4 = interpol[_x+8][_y+8].y<<8; // plus de precision dans la variation)

      _v = v = v1;                 // on initialise le x a interpoler

      dvdx   = (-v1+v2)>>3;        // variation en x quand y=0
      dvdy   = (-v1+v3)>>3;        // variation en y quand x=0
      dvdxdy = (v1-v2-v3+v4)>>6;   // variation de dvdx quand y varie

      for(y=_y;y<_y+8;y++) {       // on parcourt le carre de haut en bas
        for(x=_x;x<_x+8;x++) {     // de gauche a droite
          v += dvdx;               // on fait varier la valeur a interpoler
          interpol[x][y].y = v>>8; // on repasse la valeur sur 8 bits
        }
        _v   += dvdy;              // on fait varier la valeur de gauche
        v     = _v;                // on recommence au debut de la ligne
        dvdx += dvdxdy;            // on fait varier dvdx
      }

    }

}