/* Auteur: Nicolas JANEY                   */
/* nicolas.janey@univ-fcomte.fr            */
/* Avril 2004                              */
/* Fonction de calcul de lancer de rayons  */

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

#include "LancerDeRayonsFonctions.h"

#define DMIN 100000.0
#define EPSILON 0.0001

double sensibilite = 15.0;
float fact = 2.0F;
static int prof = 15;

void matricePlacement(objet *o,matrice m) {
  toScale(m,o->rx,o->ry,o->rz);
  matrice r;
  toRotation(r,o->a,o->ax,o->ay,o->az);
  matrice t;
  toTranslation(t,o->tx,o->ty,o->tz);
  produitMatriceMatrice(r,m,m);
  produitMatriceMatrice(t,m,m);
}

void matricePlacementInterm(objet *o,matrice m) {
  toRotation(m,o->a,o->ax,o->ay,o->az);
  matrice t;
  toTranslation(t,o->tx,o->ty,o->tz);
  produitMatriceMatrice(t,m,m);
}

void inverseMatricePlacement(objet *o,matrice m) {
  matrice s;
  toScale(s,1.0/o->rx,1.0/o->ry,1.0/o->rz);
  matrice r;
  toRotation(r,-o->a,o->ax,o->ay,o->az);
  toTranslation(m,-o->tx,-o->ty,-o->tz);
  produitMatriceMatrice(r,m,m);
  produitMatriceMatrice(s,m,m);
}

double intersectionCube(scene *s,objet *o,rayon *r) {
  o->ori[0] = r->o[0];
  o->ori[1] = r->o[1];
  o->ori[2] = r->o[2];
  o->ori[3] = 1.0;
  produitMatriceVecteur(o->mpi,o->ori,o->ori);
  o->dir[0] = r->d[0];
  o->dir[1] = r->d[1];
  o->dir[2] = r->d[2];
  o->dir[3] = 0.0;
  produitMatriceVecteur(o->mpi,o->dir,o->dir);
  int inter = 0;
  int cpt = 0;
  double dis = DMIN;
  { double d = (0.5-o->ori[2])/o->dir[2];
    double x = o->ori[0] + d*o->dir[0];
    double y = o->ori[1] + d*o->dir[1];
    if ( ( x >= -0.5 ) && ( x <= 0.5 ) && ( y >= -0.5 ) && ( y <= 0.5 ) ) {
      inter += 1;
      if ( d > EPSILON ) {
        o->norm[0] = o->norm[1] = o->norm[3] = 0.0;
        o->norm[2] = 1.0;
        dis = d;
        cpt++; } } }
  { double d = (-0.5-o->ori[2])/o->dir[2];
    double x = o->ori[0] + d*o->dir[0];
    double y = o->ori[1] + d*o->dir[1];
    if ( ( x >= -0.5 ) && ( x <= 0.5 ) && ( y >= -0.5 ) && ( y <= 0.5 ) ) {
      inter += 1;
      if ( ( d > EPSILON ) && ( d < dis ) ) {
        o->norm[0] = o->norm[1] = o->norm[3] = 0.0;
        o->norm[2] = -1.0;
        dis = d;
        cpt++; } } }
  if ( inter == 2 ) {
    if ( cpt == 0 )
      return(DMIN);
      else
      return(dis); }
  { double d = (-0.5-o->ori[1])/o->dir[1];
    double x = o->ori[0] + d*o->dir[0];
    double z = o->ori[2] + d*o->dir[2];
    if ( ( x >= -0.5 ) && ( x <= 0.5 ) && ( z >= -0.5 ) && ( z <= 0.5 ) ) {
      inter += 1;
      if ( ( d > EPSILON ) && ( d < dis ) ) {
        o->norm[0] = o->norm[2] = o->norm[3] = 0.0;
        o->norm[1] = -1.0;
        dis = d;
        cpt++; } } }
  if ( inter == 2 ) {
    if ( cpt == 0 )
      return(DMIN);
      else
      return(dis); }
  { double d = (0.5-o->ori[1])/o->dir[1];
    double x = o->ori[0] + d*o->dir[0];
    double z = o->ori[2] + d*o->dir[2];
    if ( ( x >= -0.5 ) && ( x <= 0.5 ) && ( z >= -0.5 ) && ( z <= 0.5 ) ) {
      inter += 1;
      if ( ( d > EPSILON ) && ( d < dis ) ) {
        o->norm[0] = o->norm[2] = o->norm[3] = 0.0;
        o->norm[1] = 1.0;
        dis = d;
        cpt++; } } }
  if ( inter == 2 ) {
    if ( cpt == 0 )
      return(DMIN);
      else
      return(dis); }
  { double d = (-0.5-o->ori[0])/o->dir[0];
    double y = o->ori[1] + d*o->dir[1];
    double z = o->ori[2] + d*o->dir[2];
    if ( ( y >= -0.5 ) && ( y <= 0.5 ) && ( z >= -0.5 ) && ( z <= 0.5 ) ) {
      inter += 1;
      if ( ( d > EPSILON ) && ( d < dis ) ) {
        o->norm[1] = o->norm[2] = o->norm[3] = 0.0;
        o->norm[0] = -1.0;
        dis = d;
        cpt++; } } }
  if ( inter == 2 ) {
    if ( cpt == 0 )
      return(DMIN);
      else
      return(dis); }
  { double d = (0.5-o->ori[0])/o->dir[0];
    double y = o->ori[1] + d*o->dir[1];
    double z = o->ori[2] + d*o->dir[2];
    if ( ( y >= -0.5 ) && ( y <= 0.5 ) && ( z >= -0.5 ) && ( z <= 0.5 ) ) {
      inter += 1;
      if ( ( d > EPSILON ) && ( d < dis ) ) {
        o->norm[1] = o->norm[2] = o->norm[3] = 0.0;
        o->norm[0] = 1.0;
        dis = d;
        cpt++; } } }
  if ( inter == 2 )
    if ( cpt == 0 )
      return(DMIN);
      else
      return(dis);
  return(DMIN);
}

double intersectionSphere(scene *s,objet *o,rayon *r) {
  o->ori[0] = r->o[0];
  o->ori[1] = r->o[1];
  o->ori[2] = r->o[2];
  o->ori[3] = 1.0;
  produitMatriceVecteur(o->mpi,o->ori,o->ori);
  o->dir[0] = r->d[0];
  o->dir[1] = r->d[1];
  o->dir[2] = r->d[2];
  o->dir[3] = 0.0;
  produitMatriceVecteur(o->mpi,o->dir,o->dir);
  double a = o->dir[0]*o->dir[0]+o->dir[1]*o->dir[1]+o->dir[2]*o->dir[2];
  double b = 2*(o->dir[0]*o->ori[0]+o->dir[1]*o->ori[1]+o->dir[2]*o->ori[2]);
  double c = o->ori[0]*o->ori[0]+o->ori[1]*o->ori[1]+o->ori[2]*o->ori[2] - 1.0;
  double delta = b*b - 4*a*c;
  if ( delta == 0.0 ) {
    double d = -b/2/a;
    if ( d > EPSILON ) {
      o->norm[0] = (o->dir[0]*d + o->ori[0])*o->ry*o->rz;
      o->norm[1] = (o->dir[1]*d + o->ori[1])*o->rx*o->rz;
      o->norm[2] = (o->dir[2]*d + o->ori[2])*o->rx*o->ry;
      o->norm[3] = 0.0;
      return(d); }
    else
    return(DMIN); }
  if ( delta > 0.0 ) {
    double d1 = (-b-sqrt(delta))/2/a;
    if ( d1 > EPSILON ) {
      o->norm[0] = (o->dir[0]*d1 + o->ori[0])*o->ry*o->rz;
      o->norm[1] = (o->dir[1]*d1 + o->ori[1])*o->rx*o->rz;
      o->norm[2] = (o->dir[2]*d1 + o->ori[2])*o->rx*o->ry;
      o->norm[3] = 0.0;
      return(d1); }
      else {
      double d2 = (-b+sqrt(delta))/2/a;
      if ( d2 > EPSILON ) {
        o->norm[0] = (o->dir[0]*d2 + o->ori[0])*o->ry*o->rz;
        o->norm[1] = (o->dir[1]*d2 + o->ori[1])*o->rx*o->rz;
        o->norm[2] = (o->dir[2]*d2 + o->ori[2])*o->rx*o->ry;
        o->norm[3] = 0.0;
        return(d2); } } }
  return(DMIN);
}

int testIntersectionCube(scene *s,objet *o,rayon *r) {
  vecteur ori;
  vecteur dir;
  ori[0] = r->o[0];
  ori[1] = r->o[1];
  ori[2] = r->o[2];
  ori[3] = 1.0;
  produitMatriceVecteur(o->mpi,ori,ori);
  dir[0] = r->d[0];
  dir[1] = r->d[1];
  dir[2] = r->d[2];
  dir[3] = 0.0;
  produitMatriceVecteur(o->mpi,dir,dir);
  { double d = (0.5-ori[2])/dir[2];
    double x = ori[0] + d*dir[0];
    double y = ori[1] + d*dir[1];
    if ( ( x >= -0.5 ) && ( x <= 0.5 ) && ( y >= -0.5 ) && ( y <= 0.5 ) && ( d > EPSILON ) ) {
      return(1); } }
  { double d = (-0.5-ori[2])/dir[2];
    double x = ori[0] + d*dir[0];
    double y = ori[1] + d*dir[1];
    if ( ( x >= -0.5 ) && ( x <= 0.5 ) && ( y >= -0.5 ) && ( y <= 0.5 ) && ( d > EPSILON ) ) {
      return(1); } }
  { double d = (-0.5-ori[1])/dir[1];
    double x = ori[0] + d*dir[0];
    double z = ori[2] + d*dir[2];
    if ( ( x >= -0.5 ) && ( x <= 0.5 ) && ( z >= -0.5 ) && ( z <= 0.5 ) && ( d > EPSILON ) ) {
      return(1); } }
  { double d = (0.5-ori[1])/dir[1];
    double x = ori[0] + d*dir[0];
    double z = ori[2] + d*dir[2];
    if ( ( x >= -0.5 ) && ( x <= 0.5 ) && ( z >= -0.5 ) && ( z <= 0.5 ) && ( d > EPSILON ) ) {
      return(1); } }
  { double d = (-0.5-ori[0])/dir[0];
    double y = ori[1] + d*dir[1];
    double z = ori[2] + d*dir[2];
    if ( ( y >= -0.5 ) && ( y <= 0.5 ) && ( z >= -0.5 ) && ( z <= 0.5 ) && ( d > EPSILON ) ) {
      return(1); } }
  { double d = (0.5-ori[0])/dir[0];
    double y = ori[1] + d*dir[1];
    double z = ori[2] + d*dir[2];
    if ( ( y >= -0.5 ) && ( y <= 0.5 ) && ( z >= -0.5 ) && ( z <= 0.5 ) && ( d > EPSILON ) ) {
      return(1); } }
  return(0);
}

int testIntersectionSphere(scene *s,objet *o,rayon *r) {
  vecteur ori;
  vecteur dir;
  ori[0] = r->o[0];
  ori[1] = r->o[1];
  ori[2] = r->o[2];
  ori[3] = 1.0;
  produitMatriceVecteur(o->mpi,ori,ori);
  dir[0] = r->d[0];
  dir[1] = r->d[1];
  dir[2] = r->d[2];
  dir[3] = 0.0;
  produitMatriceVecteur(o->mpi,dir,dir);
  double a = dir[0]*dir[0]+dir[1]*dir[1]+dir[2]*dir[2];
  double b = 2*(dir[0]*ori[0]+dir[1]*ori[1]+dir[2]*ori[2]);
  double c = ori[0]*ori[0]+ori[1]*ori[1]+ori[2]*ori[2] - 1.0;
  double delta = b*b - 4*a*c;
  if ( delta == 0.0 ) {
    double d = -b/2/a;
    if ( d > EPSILON ) {
      return(1); } }
  if ( delta > 0.0 ) {
    double d1 = (-b-sqrt(delta))/2/a;
    if ( d1 > EPSILON ) {
      return(1); }
      else {
      double d2 = (-b+sqrt(delta))/2/a;
      if ( d2 > EPSILON ) {
        return(1); } } }
  return(0);
}

int rayonTransmis(objet *obj,rayon *r,int interieur) {
  double indice = ( interieur ) ? obj->m->indice : 1.0/obj->m->indice;
  vecteur i = { -r->d[0],-r->d[1],-r->d[2],0.0 } ;
  vecteur norm = { obj->norm[0],obj->norm[1],obj->norm[2],0.0 } ;
  if ( interieur ) {
    norm[0] = -norm[0];
    norm[1] = -norm[1];
    norm[2] = -norm[2]; }
  double scal = produitScalaire(norm,i);
  double n2 = indice*indice;
  double val = 1.0F-n2*(1.0F-scal*scal);
  if ( val < 0.0F )
    return(0);
    else {
    val = indice*scal-sqrt(val);
    obj->rt.d[0] = val*norm[0]-indice*i[0];
    obj->rt.d[1] = val*norm[1]-indice*i[1];
    obj->rt.d[2] = val*norm[2]-indice*i[2];
    obj->rt.d[3] = 0.0;
    obj->rt.o[0] = obj->p[0];
    obj->rt.o[1] = obj->p[1];
    obj->rt.o[2] = obj->p[2];
    obj->rt.o[3] = 1.0;
    return(1); }
}

void rayonReflechi(objet *obj,rayon *r) {
  vecteur i = { -r->d[0],-r->d[1],-r->d[2],0.0 } ;
  double scal = produitScalaire(i,obj->norm);
  obj->rr.d[0] = 2*obj->norm[0]*scal-i[0];
  obj->rr.d[1] = 2*obj->norm[1]*scal-i[1];
  obj->rr.d[2] = 2*obj->norm[2]*scal-i[2];
  obj->rr.d[3] = 0.0;
  obj->rr.o[0] = obj->p[0];
  obj->rr.o[1] = obj->p[1];
  obj->rr.o[2] = obj->p[2];
  obj->rr.o[3] = 1.0;
}

void calculCaracteristiquesIntersection(objet *obj,double dmin,struct rayon *r,int interieur) {
  obj->p[0] = obj->ori[0] + dmin*obj->dir[0];
  obj->p[1] = obj->ori[1] + dmin*obj->dir[1];
  obj->p[2] = obj->ori[2] + dmin*obj->dir[2];
  obj->p[3] = 1.0;
  produitMatriceVecteur(obj->mp,obj->p,obj->p);
  produitMatriceVecteur(obj->mpp,obj->norm,obj->norm);
  normalise(obj->norm);
  if ( ( obj->m->r ) && ( !interieur ) )
    rayonReflechi(obj,r);
  if ( obj->m->t ) {
    obj->vrt = rayonTransmis(obj,r,interieur); }
    else
    obj->vrt = 0;
}

void diffusion(objet *obj,lumiere *lum,couleur *cd) {
  double dist = distance(lum->pos,obj->p);
  vecteur i;
  calculVecteurNorme(obj->p,lum->pos,i) ;
  double fact = produitScalaire(i,obj->norm);
  if ( fact > 0 ) {
    fact = fact/dist/dist*sensibilite;
    cd->r =(float) (fact*lum->coul->r*lum->e);
    cd->v =(float) (fact*lum->coul->v*lum->e);
    cd->b =(float) (fact*lum->coul->b*lum->e); }
    else {
    cd->r = cd->v = cd->b = 0.0F; }
}

int testRayonOmbre(scene *sc,objet *obj,lumiere *lum) {
  rayon r;
  calculVecteur(obj->p,lum->pos,r.d) ;
  r.o[0] = obj->p[0];
  r.o[1] = obj->p[1];
  r.o[2] = obj->p[2];
  r.o[3] = 1.0;
  r.d[3] = 0.0;
  for ( int i = 0 ; i < sc->nobj ; i++ ) {
    switch ( sc->obj[i].type ) {
      case cube      : if ( testIntersectionCube(sc,&sc->obj[i],&r) ) {
                         return(1); }
                       break;
      case sphere    : if ( testIntersectionSphere(sc,&sc->obj[i],&r) ) {
                         return(1); }
                       break; } }
  return(0);
}

void lancerDeRayons(scene *sc,rayon *r,couleur *c,int prof,int interieur) {
  double dmin = DMIN;
  int obj = -1;
  c->r = c->v = c->b = 0.0F;
  for ( int i = 0 ; i < sc->nobj ; i++ ) {
    double d;
    switch ( sc->obj[i].type ) {
      case cube      : d = intersectionCube(sc,&sc->obj[i],r);
                       break;
      case sphere    : d = intersectionSphere(sc,&sc->obj[i],r);
                       break; }
    if ( ( d > EPSILON ) && ( d < dmin ) ) {
      dmin = d;
      obj = i; } }
   if ( obj != -1 ) {
     calculCaracteristiquesIntersection(&sc->obj[obj],dmin,r,interieur);
     objet aux = sc->obj[obj];
     for ( int l = 0 ; l < sc->nlum ; l++ ) {
       if ( !testRayonOmbre(sc,&sc->obj[obj],&sc->lum[l]) ) {
         couleur cd;
         diffusion(&sc->obj[obj],&sc->lum[l],&cd);
         c->r += cd.r*sc->obj[obj].m->diffuse->r;
         c->v += cd.v*sc->obj[obj].m->diffuse->v;
         c->b += cd.b*sc->obj[obj].m->diffuse->b;
         sc->obj[obj] = aux; } }
     if ( prof > 0 ) {
       if ( sc->obj[obj].vrt ) {
         couleur ct;
         lancerDeRayons(sc,&sc->obj[obj].rt,&ct,prof-1,!interieur);
         c->r += ct.r*sc->obj[obj].m->transparence->r;
         c->v += ct.v*sc->obj[obj].m->transparence->v;
         c->b += ct.b*sc->obj[obj].m->transparence->b;
         sc->obj[obj] = aux; }
       if ( ( sc->obj[obj].m->r ) && ( !interieur ) ) {
         couleur cr;
         lancerDeRayons(sc,&sc->obj[obj].rr,&cr,prof-1,0);
         if ( sc->obj[obj].vrt ) {
           c->r += cr.r*sc->obj[obj].m->specular->r;
           c->v += cr.v*sc->obj[obj].m->specular->v;
           c->b += cr.b*sc->obj[obj].m->specular->b; }
           else {
           c->r += cr.r*(sc->obj[obj].m->specular->r+sc->obj[obj].m->transparence->r);
           c->v += cr.v*(sc->obj[obj].m->specular->v+sc->obj[obj].m->transparence->v);
           c->b += cr.b*(sc->obj[obj].m->specular->b+sc->obj[obj].m->transparence->b); } } } }
}

void initialisationsCalculLancerDeRayon(scene *sc) {
  for ( int o = 0 ; o < sc->nobj ; o++ ) {
    inverseMatricePlacement(&sc->obj[o],sc->obj[o].mpi);
    matricePlacementInterm(&sc->obj[o],sc->obj[o].mpp);
    matricePlacement(&sc->obj[o],sc->obj[o].mp); }
}

unsigned char *initialisationsCalculImageLancerDeRayon(scene *sc,int width2,int height2) {
  initialisationsCalculLancerDeRayon(sc);
  int w = width2;
  int h = height2;
  int t = w*h*3 ;
  return((unsigned char *) calloc(t,sizeof(unsigned char)));
}

void calculImageLancerDeRayon(scene *sc,int j,int w,int h,unsigned char *im) {
  double d = tan(sc->vang/360.0F*3.14159F)*sc->cmin/h*2.0F;
  for ( int i = 0 ; i < w ; i++ ) {
    int p = i*3;
    rayon r;
    r.o[0] = 0.0F;
    r.o[1] = 0.0F;
    r.o[2] = 0.0F;
    r.o[3] = 1.0F;
    r.d[0] = (i-w/2)*d+d/2.0F;
    r.d[1] = (j-h/2)*d+d/2.0F;
    r.d[2] = -sc->cmin;
    r.d[3] = 1.0F;
    normalise(r.d);
    couleur c;
    lancerDeRayons(sc,&r,&c,prof,0);
    if ( c.r > 1.0F )
      c.r = 1.0F;
    if ( c.v > 1.0F )
      c.v = 1.0F;
    if ( c.b > 1.0F )
      c.b = 1.0F;
    im[p] =(unsigned char) (c.r*255.0F);
    im[p+1] =(unsigned char) (c.v*255.0F);
    im[p+2] =(unsigned char) (c.b*255.0F); }
}