#include <GL/glut.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "types.h"
#include "globj.h"
#include "vecmath.h"
#include "sphere.h"
#include "laminate.h"
#include "glcam.h"
#include "gloutils.h"

extern int selectionx, selectiony ;
void drawsphere    (globj_t*, glcam_t* camera) ;
void maptosphere (globj_t*, vector*) ;
void sphere_init () ;
globj_t* sphere_factory () ;
void destructsphere (globj_t*);
int recursion = 8 ;

globj_class_t globj_class_sphere = 
{ sphere_init, sphere_factory, "sphere", NULL } ;

globj_t* sphere_factory ()
{
  globj_t* globj = malloc (sizeof(globj_t)) ;
  globj->drawfunc = &drawsphere ;
  globj->destructor = &destructsphere ;
  globj->maptosurface = &maptosphere ;
  globj->name = "sphere" ;
  globj->data = sphere_new () ;
  return globj ;
}

void maptosphere (globj_t* globj, vector* point)
{
  sphere_maptosurface (globj->data, point) ;
}

void drawsphere (globj_t* globj, glcam_t* camera)
{
  //sphere_drawfast (globj->data) ;
  sphere_drawopt(globj->data, camera, recursion, 100);
}

void destructsphere (globj_t* globj)
{
  sphere_delete (globj->data) ;
  free (globj) ;
}

/* partie concernant glhand */
int detailCallback (void* data, int eventtype, int eventdata, int x, int y)
{
  sphere_t* sphere = *((sphere_t**)data) ;
  sphere_detail (sphere) ;
  return 1;
}

int laminateCallback (void* data, int eventtype, int eventdata, int x, int y)
{
  sphere_t* sphere = *((sphere_t**)data) ;
  sphere_adjust (sphere) ;
  return 1 ;
}

int recursionCallback (void* data, int eventtype, int eventdata, int x, int y)
{
  int* reclevel = data ;
  switch (eventdata) {
  case 'R': (*reclevel)++; return 1;
  case 'r': (*reclevel)--; return 1;
  default: return 0;
  }
}

action_t action_sphere_detail = 
{ "detail", detailCallback, &zone_sphere.data } ;
action_t action_sphere_laminate = 
{ "laminate", laminateCallback, &zone_sphere.data } ;
action_t action_sphere_recursion =
{ "recursion", recursionCallback, &recursion } ;

mapping_t mappings_sphere [] = {
  { EVENT_KEY, 'd', &action_sphere_detail },
  { EVENT_KEY, 'l', &action_sphere_laminate },
  { EVENT_KEY, 'r', &action_sphere_recursion },
  { EVENT_KEY, 'R', &action_sphere_recursion }
};

zone_t zone_sphere = { 0, NULL, sizeof (mappings_sphere)/sizeof(mapping_t),
		       mappings_sphere, NULL } ;
/* fin interface glhand */

void sphere_init (void) {}

#define X .525731112119133606 
#define Z .850650808352039932

/*
static vector icosa_vertices[12] = {    
  {-X, 0.0, Z}, {X, 0.0, Z}, {-X, 0.0, -Z}, {X, 0.0, -Z},    
  {0.0, Z, X}, {0.0, Z, -X}, {0.0, -Z, X}, {0.0, -Z, -X},    
  {Z, X, 0.0}, {-Z, X, 0.0}, {Z, -X, 0.0}, {-Z, -X, 0.0} 
};

static int icosa_triangles[20][3] = { 
  {0,1,4}, {0,4,9}, {9,4,5}, {4,8,5}, {4,1,8},     
  {8,1,10}, {8,10,3}, {5,8,3}, {5,3,2}, {2,3,7},    
  {7,3,10}, {7,10,6}, {7,6,11}, {11,6,0}, {0,6,1}, 
  {6,10,1}, {9,11,0}, {9,2,11}, {9,5,2}, {7,11,2} 
};

static int icosa_links[12][5] = {
  {0,1,13,14,16}, {0,4,5,6,7}, {8,9,17,18,19}, {6,7,8,9,10},
  {0,1,2,3,4}, {3,4,7,8,18}, {11,12,13,14,15}, {9,10,11,12,19},
  {3,4,5,6,7}, {1,2,16,17,18}, {5,6,10,11,15}, {12,13,16,17,19}
};
*/

static vector octa_vertices[6] = {
  {0.0, 0.0, 0.5}, {-0.5, 0.0, 0.0}, {0.0, 0.0, -0.5}, {0.5, 0.0, 0.0}, {0.0, 0.5, 0.0}, {0.0, -0.5, 0.0}
};

static int octa_triangles[8][3] = {
  {0,4,1}, {1,4,2}, {2,4,3}, {3,4,0},
  {0,5,3}, {3,5,2}, {2,5,1}, {1,5,0}
};

static int octa_links[6][5] = {
  {0,3,4,7}, {0,1,6,7}, {1,2,5,6}, {2,3,4,5}, {0,1,2,3}, {4,5,6,7}
};

sphere_t* sphere_new ()
{
  int i,j;
  sphere_t* sphere    = malloc (sizeof(sphere_t)) ;
  sphere->vertices    = malloc (6*sizeof(vertex_t)) ;
  sphere->triangles   = malloc (8*sizeof(triangle_t)) ;
  sphere->n_vertices  = 6 ;
  sphere->n_triangles = 8 ;
  sphere->n_internals = 0 ;
  sphere->n_lines     = 12 ;
  sphere->n_roots     = 8 ;
  sphere->depth       = 0 ;
  for (i=0; i<6; i++) {
    sphere->vertices[i].position = octa_vertices[i] ;
    vector_assign(&sphere->vertices[i].normal, 0.0, 0.0, 0.0);
    for (j=0; j<5; j++) 
      sphere->vertices[i].triangles[j] = octa_links[i][j] ;
  }
  for (i=0; i<8; i++) {
    for (j=0; j<3; j++) 
      sphere->triangles[i].vertices[j] = octa_triangles[i][j] ;
    for (j=0; j<4; j++)
      sphere->triangles[i].triangles[j] = -1 ;
    vector_computenormal 
      (&sphere->triangles[i].normal,
       sphere->vertices[sphere->triangles[i].vertices[0]].position,
       sphere->vertices[sphere->triangles[i].vertices[1]].position,
       sphere->vertices[sphere->triangles[i].vertices[2]].position) ;
    sphere->triangles[i].root = -1 ;
  }
  
  for (i=0; i<6; i++)
    vertex_setcolor (sphere->vertices+i, 0.9, 0.7, 0.3, 1.0) ;

  for (i=0; i<6; i++)
    sphere_computevertexnormal(&sphere->vertices[i], sphere);
  
  zone_sphere.data = sphere ;
  return sphere ;
}

void sphere_delete (sphere_t* sphere)
{
  free (sphere->vertices) ;
  free (sphere->triangles) ;
  free (sphere) ;
}

void sphere_drawfull (sphere_t* sphere, int drawinternals)
{
  int i,j ;
  if (drawinternals) i=0 ; else i=sphere->n_internals;
  for (; i<sphere->n_triangles; i++) {
    triangle_t* triangle = &sphere->triangles[i] ;
    glNormal3fv ((GLfloat*)&triangle->normal) ;
    glBegin (GL_TRIANGLES) ;
    for (j=0; j<3; j++) {
      glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, 
		    sphere->vertices[triangle->vertices[i]].color) ;
      glVertex3fv ((GLfloat*)&sphere->vertices[triangle->vertices[i]].position);
    }
    glEnd() ; 
  }
}

/* draws only LEAF triangles */
void sphere_drawfast (sphere_t* sphere)
{
  sphere_drawfull (sphere, 0) ;
}

/* draws EVERY triangle */
void sphere_drawslow (sphere_t* sphere)
{
  sphere_drawfull (sphere, 1) ;
}


/* draws only visible triangles */
void sphere_drawopt (sphere_t* sphere, glcam_t* camera, int depth, int minsize)
{
  opengldata_t gldata ;
  sphere_drawtriangle_f drawtriangle = NULL;
  int i,count=0 ;
  matrix inverse;
  GLint shademodel, rendermode ;

  glouGet (&gldata) ;
  matrix_inverse (&inverse, gldata.modelviewmatrix) ;

  glDisableClientState (GL_VERTEX_ARRAY) ;  
  glDisableClientState (GL_NORMAL_ARRAY) ;
  glDisableClientState (GL_COLOR_ARRAY) ;
  glDisable (GL_COLOR_MATERIAL) ;

  glEnableClientState (GL_VERTEX_ARRAY) ;
  glVertexPointer (3, GL_COORD, sizeof(vertex_t), sphere->vertices) ;

  /* examine openGL state, and determine which triangle rasterizer to call */
  
  glGetIntegerv (GL_RENDER_MODE, &rendermode) ;
  switch (rendermode) {
  case GL_RENDER:
    glGetIntegerv (GL_SHADE_MODEL, &shademodel) ;
    if (shademodel == GL_SMOOTH) {
      glEnableClientState (GL_COLOR_ARRAY) ;
      glColorPointer (4,GL_FLOAT,sizeof(vertex_t), &(sphere->vertices->color));
      if (glIsEnabled (GL_LIGHTING)) {
	drawtriangle = sphere_drawtriangle_gl_smoothlight ;
	glEnableClientState (GL_NORMAL_ARRAY) ;
	glNormalPointer (GL_COORD, sizeof(vertex_t),&(sphere->vertices->normal));
	glColorMaterial (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE) ;
	glEnable (GL_COLOR_MATERIAL) ;
      } else { /* GL_SMOOTH, lighting disabled */
	drawtriangle = sphere_drawtriangle_gl_smoothcolor ;
      }
    }
    else /* GL_FLAT */
      {
	if (glIsEnabled (GL_LIGHTING))
	  drawtriangle = sphere_drawtriangle_gl_flatlight ;
	else
	  drawtriangle = sphere_drawtriangle_gl_flatcolor ;
      }
    break;
  case GL_SELECT:
    glMatrixMode (GL_PROJECTION) ;
    glLoadIdentity () ;
    gluPickMatrix (selectionx, gldata.viewport[3]-selectiony, 
		   4, 4, gldata.viewport) ;
    glMultMatrixd (gldata.projection) ;
    glMatrixMode (GL_MODELVIEW) ;
    drawtriangle = sphere_drawtriangle_gl_selection ;
    break;
  default:
    fprintf (stderr, "Unknown RENDERMODE!\n") ;
  }
  
  glInitNames(); /* pour la selection */
  glPushName(-1) ;
  glBegin (GL_TRIANGLES) ;
  for(i=0; i<8; i++) {
    sphere_drawtriangle(sphere,sphere->triangles+i, &gldata, camera, 
			depth, minsize, &count, drawtriangle);
  }
  glEnd () ;
  glLoadName(-1) ;
  glDisableClientState (GL_VERTEX_ARRAY) ;  
  glDisableClientState (GL_NORMAL_ARRAY) ;
  glDisableClientState (GL_COLOR_ARRAY) ;
  glDisable (GL_COLOR_MATERIAL) ;
  
  {
    char text [] = "%d faces" ;
    int i,len=strlen(text)+10 ;
    char buffer [len] ;
    snprintf (buffer, len, text, count) ;
    glPushAttrib (GL_LIGHTING_BIT) ;
    glDisable (GL_LIGHTING) ;
    glColor3f (1,1,1) ;
    glouAbsoluteRasterPos3f (1, -24, 0.1) ;
    for (i = 0; i < len && buffer[i]; i++) {
      glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_24, buffer[i]);
    }
    glPopAttrib () ;
  }

}

void sphere_drawtriangle_gl_highlight (sphere_t* s, triangle_t* t)
{
  int i ;
  for (i=0; i<3; i++)
    glVertex3fv ((GLfloat*)&s->vertices[t->vertices[i]].position) ;
}

void sphere_drawtriangle_gl_selection (sphere_t* s, triangle_t* t)
{
  int i ;
  glEnd() ;
  glLoadName (sphere_findtriangleindex (s,t)) ;
  glBegin (GL_TRIANGLES) ;
  for (i=0; i<3; i++)
    glArrayElement (t->vertices[i]) ;
}

void sphere_drawtriangle_gl_flatcolor (sphere_t* s, triangle_t* t)
{
  int i,j ;
  GLfloat color[4] ;
  for (i=0; i<4; i++) {
    color[i] = 0;
    for (j=0; j<3; j++)
      color[i] += s->vertices[t->vertices[j]].color[i] ;
    color[i] /= 3 ;
  }
  glColor4fv (color) ;
  for (i=0; i<3; i++)
    glArrayElement (t->vertices[i]) ;
}

void sphere_drawtriangle_gl_smoothcolor (sphere_t* s, triangle_t* t)
{
  int i ;
  for (i=0; i<3; i++)
    glArrayElement (t->vertices[i]) ;
}

void sphere_drawtriangle_gl_smoothlight (sphere_t* s, triangle_t* t)
{
  int i ;
  for (i=0; i<3; i++)
    glArrayElement (t->vertices[i]) ;
}

void sphere_drawtriangle_gl_flatlight (sphere_t* s, triangle_t* t)
{
  int i,j ;
  GLfloat color[4] ;
  //vector normal ;
  for (i=0; i<4; i++) {
    color[i] = 0;
    for (j=0; j<3; j++)
      color[i] += s->vertices[t->vertices[j]].color[i] ;
    color[i] /= 3 ;
  }
  /*
    vector_computenormal (&normal, 
    s->vertices[t->vertices[0]].normal,
    s->vertices[t->vertices[2]].normal,
    s->vertices[t->vertices[1]].normal) ;
    glNormal3fv ((GLfloat*)&normal) ;
  */
  glNormal3fv ((GLfloat*)&t->normal) ;
  glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, color) ;
  for (i=0; i<3; i++)
    glArrayElement (t->vertices[i]) ;
}

void sphere_drawsubtriangles (sphere_t* sphere, triangle_t* triangle,
			      opengldata_t* gldata, glcam_t* cam, 
			      int depth, int minsize, int* count,
			      sphere_drawtriangle_f drawtriangle)
{
  int i ;
  for (i=0; i<4; i++)
    sphere_drawtriangle (sphere, sphere->triangles+triangle->triangles[i], 
			 gldata, cam, depth, minsize, count, drawtriangle) ;
}

void sphere_drawtriangle(sphere_t* sphere, triangle_t* triangle,
			 opengldata_t* gldata, glcam_t* camera, 
			 int depth, int minsize, int* count,
			 sphere_drawtriangle_f drawtriangle)
{
  int size,i,visiblecount=0 ;
  GLdouble x[3],y[3],z[3] ;

  if (camera) {
    /* essayons donc d'optimiser en ne dessinant que ce qu'on voit ! */
    vector eye = camera->eye;
    vector_applymatrix (&eye, gldata->modelviewinverse, 1.0) ;
    for (i=0; i<3; i++) { /* testons les 3 normales des sommets de la face */
      vector dir = sphere->vertices[triangle->vertices[i]].position;
      vector n;
      vector_sub(&dir, eye);
      n = sphere->vertices[triangle->vertices[i]].normal;
      if (vector_scalar(dir, n) < 0.1) visiblecount++ ;
    }
    if (visiblecount==0) return ; /* aucune normale "visible", on arrete */
    if (visiblecount==3) camera=NULL; /* toutes normales "visibles", tout aff*/
  }

  /* check whether triangle is on screen or not */
  visiblecount=0 ;
  for (i=0; i<3; i++) {
    gluProject (sphere->vertices[triangle->vertices[i]].position.x,
		sphere->vertices[triangle->vertices[i]].position.y,
		sphere->vertices[triangle->vertices[i]].position.z,
		gldata->modelview, gldata->projection, gldata->viewport,
		&x[i], &y[i], &z[i]) ;
    if (x[i]<gldata->viewport[0]) continue ;
    if (x[i]>gldata->viewport[0]+gldata->viewport[2]) continue ;
    if (y[i]<gldata->viewport[1]) continue ;
    if (y[i]>gldata->viewport[1]+gldata->viewport[3]) continue ;
    visiblecount++ ;
  }
  //if (visiblecount==0) return ; /* triangle completely offscreen */
  /* FIXME: this does not account for triangles with all vertices offscreen,
     but crossing the screen anyway! */
  /* FIXME: cache visiblecount results */

  //size = triangle_computearea (x,y) ;
  //if (size<minsize) depth=0 ;

  /* arrive ici, soit on avait cam=NULL (pas d'optimisation),
     soit on a determine qu'il fallait afficher la face et ses filles */
  if (depth && triangle->triangles[0] != -1) { /* faut-il recursiver ? */
    sphere_drawsubtriangles (sphere, triangle, gldata, camera, 
			     depth-1, minsize, count, drawtriangle) ;
    return ;
  }

  /* quand on arrive la, c'est qu'il faut dessiner la face */
  (*count)++;
  drawtriangle (sphere, triangle) ;
}

typedef struct newvertex_s {
  int neighbour, middle ;
} newvertex_t ;
newvertex_t newvertex = { -1, -1 } ;

void vertex_interpolate (vertex_t* vnew, vertex_t* v1, vertex_t* v2)
{
  /* we interpolate vnew as if it were located on a 0-centered spheroid,
     between v1 and v2 .  */
  coord h ; /* height of the vertex exactly between v1 & v2 */
  coord l ; /* desired height of vnew */
  int i;

  vnew->position = v1->position ;
  vector_add (&vnew->position, v2->position) ;
  h = vector_norm (&vnew->position) ;
  l = 0.5 * (vector_norm (&v1->position) + vector_norm (&v2->position)) ;
  vector_scale (&vnew->position, l/h) ;

  /* we also interpolate colors */
  for (i=0; i<4; i++) vnew->color[i] = 0.5 * (v1->color[i]+v2->color[i]) ;
}

/* lookup vertex between (startvertex,endvertex), eventually adding it
   to the cache, and returning its index */
int vertexcache_addline (newvertex_t** vertexcache,
			  int* firstfreevertex,
			  int startvertex, int endvertex)
{
  int i,start,end ;

  /* middle of (start,end) when start==end ? guess what ... */
  if (startvertex == endvertex) return startvertex ;

  /* vcache entries are always stored as (start,end) pairs, with start<end */
  if (startvertex > endvertex) 
    { start=endvertex; end=startvertex; }
  else
    { start=startvertex; end=endvertex; }
  
  /* lookup vertices pairs starting with 'start' */
  for (i=0; (i<6) && (vertexcache[start][i].neighbour != end); i++)
    if (vertexcache[start][i].neighbour == -1) { /* free slot found */
      vertexcache[start][i].neighbour = end ; /* insert vertex */
      vertexcache[start][i].middle = (*firstfreevertex)++ ;
      return *firstfreevertex - 1 ; /* return what we just inserted */
    }
  if (vertexcache[start][i].neighbour != end) /* not found, and no free slot */
    fprintf (stderr, "ERROR: can't add start=%d end=%d in vcache\n",start,end);
  return vertexcache[start][i].middle ; /* cache hit ! */
}

void vertexcache_addtriangle (newvertex_t** vertexcache, 
			      int* firstfreevertex,
			      triangle_t* triangles, int father, int son,
			      int v1a,int v1b,int v2a,int v2b,int v3a,int v3b)
{
  triangle_t* f = triangles+father ;
  triangle_t* s = triangles+son ;
  int i ;
  int v1ai = f->vertices[v1a] ;
  int v1bi = f->vertices[v1b] ;
  int v2ai = f->vertices[v2a] ;
  int v2bi = f->vertices[v2b] ;
  int v3ai = f->vertices[v3a] ;
  int v3bi = f->vertices[v3b] ;
  int v1 = vertexcache_addline (vertexcache, firstfreevertex, v1ai, v1bi) ;
  int v2 = vertexcache_addline (vertexcache, firstfreevertex, v2ai, v2bi) ;
  int v3 = vertexcache_addline (vertexcache, firstfreevertex, v3ai, v3bi) ;

  s->vertices[0] = v1 ;
  s->vertices[1] = v2 ;
  s->vertices[2] = v3 ;
  s->normal = f->normal ; /* will be recomputed later */

  for (i=0; i<4; i++) s->triangles[i] = -1 ;
  s->root = father ;
  for (i=0; i<4; i++)
    if (f->triangles[i] == -1) {
      f->triangles[i] = son ;
      return ;
    }
  fprintf (stderr, "ERROR: father=%d is full and can't add son=%d\n",
	   father, son) ;
}

void vertex_addtriangleref (vertex_t* vertex, int triangle)
{
  int i ;
  for (i=0; i<6; i++)
    if (vertex->triangles[i] == -1) {
      vertex->triangles[i] = triangle ;
      return ;
    }
  fprintf (stderr, "ERROR: vertex full, can't link triangle=%d\n", triangle) ;
}


/* refines sphere */
/* we use a vertex cache to store, for each pair (v1,v2), 
   the new vertex between them. NB : v1 < v2 (regarding indices) */
void sphere_detail (sphere_t* sphere)
{
  int i,j ;
  int n_newvertices = sphere->n_vertices + sphere->n_lines ;
  int n_newtriangles = sphere->n_triangles 
    + 4 * (sphere->n_triangles - sphere->n_internals) ;
  int n_newlines = (sphere ->n_lines * 2) + (3 * (sphere -> n_triangles - sphere -> n_internals)) ;
  newvertex_t** vertexcache ;
  vertex_t* newvertices ;
  triangle_t* newtriangles ;
  int firstfreevertex = sphere->n_vertices ;

  // 2KILL
  printf ("L(%d) = %d ; V(%d) = %d ; F(%d) = %d\n", sphere->depth+1, n_newlines, sphere->depth+1, n_newvertices, sphere->depth+1, n_newtriangles) ;

  /* initialize vertex cache */
  vertexcache = malloc (sphere->n_vertices * sizeof(newvertex_t*)) ;
  for (i=0; i<sphere->n_vertices; i++) {
    vertexcache[i] = malloc (6*sizeof(newvertex_t)) ;
    for (j=0; j<6; j++)
      vertexcache[i][j] = newvertex ;
  }

  /* create new vertex&triangle arrays */
  newvertices = malloc (n_newvertices*sizeof(vertex_t)) ;
  newtriangles = malloc (n_newtriangles*sizeof(triangle_t)) ;

  /* copy old vertices */
  for (i=0; i<sphere->n_vertices; i++)
    newvertices[i] = sphere->vertices[i] ;

  /* clear vertices triangles references */
  for (i=0; i<n_newvertices; i++)
    for (j=0; j<6; j++)
      newvertices[i].triangles[j] = -1 ;

  /* copy old triangles, unchanged */
  for (i=0; i<sphere->n_triangles; i++)
    newtriangles[i] = sphere->triangles[i] ;

  /* for each _external_ triangle, create 4 new ones */
  for (i=sphere->n_internals; i<sphere->n_triangles; i++) {
    int offset = sphere->n_triangles + 4*(i-sphere->n_internals) ;
    for (j=0; j<3; j++) /* create 3 "border" triangles */
      vertexcache_addtriangle (vertexcache, &firstfreevertex,
			       newtriangles, i, offset+j,
			       j,0, j,1, j,2) ;
    /* create "central" triangle */
    vertexcache_addtriangle (vertexcache, &firstfreevertex,
			     newtriangles, i, offset+3,
			     1,2, 2,0, 0,1) ;
  }
			 
  /* sanity check */
  if (firstfreevertex != n_newvertices)
    fprintf (stderr, "ERROR: ffv=%d != nnv=%d\n", 
	     firstfreevertex, n_newvertices) ;
  
  /* add vertex cache to new vertices */
  for (i=0; i<sphere->n_vertices; i++)
    for (j=0; (j<6) && (vertexcache[i][j].neighbour!=-1) ; j++)
      vertex_interpolate (newvertices+vertexcache[i][j].middle,
			  newvertices+vertexcache[i][j].neighbour,
			  newvertices+i) ;

  /* parse new triangles, and update vertices->triangles references */
  for (i=sphere->n_triangles; i<n_newtriangles; i++)
    for (j=0; j<3; j++)
      vertex_addtriangleref (newvertices+(newtriangles[i].vertices[j]), i) ;

  /* Recompute Normals */
  for (i=sphere->n_triangles; i<n_newtriangles; i++) {
    vector_computenormal 
      (&newtriangles[i].normal,
       newvertices[newtriangles[i].vertices[0]].position,
       newvertices[newtriangles[i].vertices[1]].position,
       newvertices[newtriangles[i].vertices[2]].position) ;
  }

  /* delete vertex cache */
  for (i=0; i<sphere->n_vertices; i++)
    free (vertexcache[i]) ;

  free (vertexcache) ;

  /* delete old triangles & vertices */
  free (sphere->triangles) ;
  free (sphere->vertices) ;

  /* update sphere with new data */
  sphere->triangles = newtriangles ;
  sphere->vertices = newvertices ;
  sphere->n_internals = sphere->n_triangles ;
  sphere->n_triangles = n_newtriangles ;
  sphere->n_vertices = n_newvertices ;
  sphere->n_lines = n_newlines ;
  sphere->depth++ ;
  /* jerome> je crois qu'il faut aussi faire le calcul pour les vertices */
  /* ca ne resoud pas le probleme mais au moins des faces apparaissent */
  /*
  for (i=0; i<sphere->n_vertices; i++) {
    sphere_computevertexnormal (sphere->vertices+i, sphere) ;
  }
  */
  sphere_computenormals(sphere);

}

/*
void sphere_computenormals(sphere_t* sphere)
{
  int i, j, k, ind;
  vertex_t* vertex;
  
  for(i=0; i<sphere->n_vertices; i++) {
    int n_faces = 4;
    
    vertex = sphere->vertices+i;
    vector_assign(&vertex->normal, 0.0, 0.0, 0.0);


    printf("%d\n", vertex->triangles[5]);
    if(vertex->triangles[5] != -1) n_faces=6;
    for (j=0 ; j<n_faces ; j++) {
      
      for(k=0; k<3; k++) {
	ind = sphere->triangles[vertex->triangles[j]].vertices[k];
	if(ind != i) {
	  vector a = vertex->position;
	  vector b = sphere->vertices[ind].position;
	  
	  vector_sub(&a, b);
	  vector_normalize(&a);

	  vector_add(&vertex->normal, a);
	}
      }
    }
    
    vector_normalize(&vertex->normal);
    
    if(vector_scalar(vertex->normal, vertex->position) < 0.0)
      vector_scale(&vertex->normal, -1);
  }
}
*/

void sphere_computevertexnormal(vertex_t* vertex, sphere_t *sphere)
{
  int i,n_vertices=4;
  vector_assign(&vertex->normal, 0.0, 0.0, 0.0);
  if(vertex->triangles[4] != -1) n_vertices=6;
  for (i=0 ; i<n_vertices ; i++)
    vector_add(&vertex->normal, 
	       sphere->triangles[vertex->triangles[i]].normal);
  
  vector_normalize (&vertex->normal);
}

coord sphere_distancetotriangle (sphere_t* sphere, int triangle, vector point)
{
  triangle_t* t = sphere->triangles+triangle ;
  vector center =      sphere->vertices[t->vertices[0]].position ;
  vector_add (&center, sphere->vertices[t->vertices[1]].position) ;
  vector_add (&center, sphere->vertices[t->vertices[2]].position) ;
  vector_scale (&center, 1.0/3.0) ;
  vector_sub (&center, point) ;
  return vector_norm (&center) ;
}

int sphere_findboundingtriangle (sphere_t* sphere, vector point)
{
  int i=0, besttriangle=-1 ;
  coord bestvalue = COORD_INFINITY ;

  /* on détermine le "toplevel" triangle au dessus duquel est le point */
  for (i=0; i<8; i++) {
    coord value = sphere_distancetotriangle (sphere, i, point) ;
    if (value<bestvalue) {
      value=bestvalue;
      besttriangle=i;
    }
  }

  /* on descend récursivement jusqu'à tomber sur une feuille */
  while (sphere->triangles[besttriangle].triangles[0] != -1) {
    int firstson = sphere->triangles[besttriangle].triangles[0] ;
    bestvalue = COORD_INFINITY ;
    for (i=0; i<3; i++) {
      coord value = sphere_distancetotriangle (sphere, firstson+i, point) ;
      if (value<bestvalue) {
	value=bestvalue;
	besttriangle = firstson+i ;
      }
    }
  }
  return besttriangle ;
}

void sphere_maptosurface (sphere_t* sphere, vector* point)
{
  int triangle = sphere_findboundingtriangle (sphere, *point) ;
  int* vertices = sphere->triangles[triangle].vertices ;
  *point =           sphere->vertices[vertices[0]].position ;
  vector_add (point, sphere->vertices[vertices[1]].position) ;
  vector_add (point, sphere->vertices[vertices[2]].position) ;
  vector_scale (point, 1.0/3.0) ;
  /* FIXME: maybe we should simulate sphere curvature. */
}

void vertex_setcolor (vertex_t* vertex, float r, float g, float b, float a) {
  vertex->color[0] = r ;
  vertex->color[1] = g ;
  vertex->color[2] = b ;
  vertex->color[3] = a ;
}

int triangle_computearea (double x[3], double y[3]) {
  double x1 = x[1]-x[0] ;
  double x2 = x[2]-x[0] ;
  double y1 = y[1]-y[0] ;
  double y2 = y[2]-y[0] ;
  return fabs(0.5*(x1*y2-x2*y1)) ;
}
 
int sphere_findtriangleindex (sphere_t* sphere, triangle_t* triangle) {
  int i ;
  if (triangle->root == -1) { /* toplevel triangle, brute force */
    for (i=0; i<sphere->n_roots; i++)
      if (triangle == sphere->triangles+i) return i;
    fprintf (stderr, "WARNING: triangle index not found!\n") ;
    return -1 ;
  }
  /* clever method : who am I ? I am a son of my father! */
  for (i=0; i<4; i++) {
    int index = sphere->triangles[triangle->root].triangles[i] ;
    if (sphere->triangles+index == triangle) return index ;
  }
  fprintf (stderr, "WARNING: triangle index not found!\n") ;
  return -1 ;
}