double calculCoordButterFly(int A, int B, int C, int D, int E, int F, int G, int H) {

  return ((TabVertices[A]/2) + TabVertices[B]/2 + TabVertices[C]/8 + TabVertices[D]/8 - TabVertices[G]/16 - TabVertices[H]/16 - TabVertices[F]/16 - TabVertices[E]/16);
}


void butterFly() {

  int oldNbPolys = nbPolys;
  nbPolys=nbPolys + 7*nbPolys;

  int oldNbPoints = nbPoints;
  nbPoints= nbPoints + 3*oldNbPolys;

  TabVertices = (float*) realloc(TabVertices, sizeof(float) * (3*nbPoints));
  TabIndices =  (int*)   realloc(TabIndices,  sizeof(int) * (3*nbPolys));
  TabNormals =  (float*) realloc(TabNormals,  sizeof(int) * (3*nbPolys));
  TabVoisins = (int*)   realloc(TabVoisins, sizeof(int)*(3*nbPolys));

  // Crée les sous triangles du centre
  int k = oldNbPolys;
  for(int i = 0; i< oldNbPolys; i++) {

    TabIndices[3*k] = oldNbPoints+3*i;
    TabIndices[3*k+1] = oldNbPoints+3*i+1;
    TabIndices[3*k+2] = oldNbPoints+3*i+2;

    std::cout<< TabIndices[3*k] <<std::endl;
    std::cout<< TabIndices[3*k+1]<<std::endl;
    std::cout<< TabIndices[3*k+2]<<std::endl<<std::endl;
    k++;
  }


  double x1, x2, x3;
  double y1, y2, y3;
  double z1, z2, z3;

  k = 2* oldNbPolys;

  for(int i = 0 ; i< oldNbPolys; i++) {

    // Pour chaque arrête du poly en question il faut faire le calcule suivant :
    // Chercher les deux sommets qui font partie des triangles qui contiennent A et B en sommets


   // premier sommet
    int A = TabIndices[3*i];
    int B = TabIndices[3*i+1];
    int C = TabIndices[3*i+2];
    std::cout<<"Point A :"<<A<<" point B : "<<B<<" point C : "<<C<<"\n";

    int D = -1;
    int G = -1;
    int H = -1;
    int F = -1;
    int E = -1;

    D = indicePointRestantTriangle(A,B, oldNbPolys, C); // Point D, symétrique de C par rapport à AB
    G = indicePointRestantTriangle(D,B, oldNbPolys, A); // Point G, symétrique de A par rapport à DB
    H = indicePointRestantTriangle(A,D, oldNbPolys, B); // Point H, symétrique de B par rapport à AD
    F = indicePointRestantTriangle(A,C, oldNbPolys, B); // Point F, symétrique de B par rapport à AC
    E = indicePointRestantTriangle(C,B, oldNbPolys, A); // Point E, symétrique de A par rapport à CB

  x1 = calculCoordButterFly(3*(A) , 3*(B), 3*C, 3*D, 3*E, 3*F, 3*G, 3*H);
  y1 = calculCoordButterFly(3*(A)+1, 3*(B)+1, 3*C+1, 3*D+1, 3*E+1, 3*F+1, 3*G+1, 3*H+1);
  z1 = calculCoordButterFly(3*(A)+2, 3*(B)+2, 3*C+2, 3*D+2, 3*E+2, 3*F+2, 3*G+2, 3*H+2);

  TabVertices[3*TabIndices[3*(oldNbPolys+i)]] =   x1;
  TabVertices[3*TabIndices[3*(oldNbPolys+i)]+1] = y1;
  TabVertices[3*TabIndices[3*(oldNbPolys+i)]+2] = z1;


  TabIndices[3*k]   = TabIndices[3*(oldNbPolys+i)]; // Ajout triangle (A,B,MILIEU)
  TabIndices[3*k+1] = A;
  TabIndices[3*k+2] = B;
  k++;

    // second sommet
    A = TabIndices[3*i];
    B = TabIndices[3*i+2];
    C = TabIndices[3*i+1];
    std::cout<<"Point A :"<<A<<" point B : "<<B<<" point C : "<<C<<"\n";

    D = -1;
    G = -1;
    H = -1;
    F = -1;
    E = -1;

    D = indicePointRestantTriangle(A,B, oldNbPolys, C); // Point D, symétrique de C par rapport à AB
    G = indicePointRestantTriangle(D,B, oldNbPolys, A); // Point G, symétrique de A par rapport à DB
    H = indicePointRestantTriangle(A,D, oldNbPolys, B); // Point H, symétrique de B par rapport à AD
    F = indicePointRestantTriangle(A,C, oldNbPolys, B); // Point F, symétrique de B par rapport à AC
    E = indicePointRestantTriangle(C,B, oldNbPolys, A); // Point E, symétrique de A par rapport à CB


  x2 = calculCoordButterFly(3*(A) , 3*(B), 3*C, 3*D, 3*E, 3*F, 3*G, 3*H);
  y2 = calculCoordButterFly(3*(A)+1, 3*(B)+1, 3*C+1, 3*D+1, 3*E+1, 3*F+1, 3*G+1, 3*H+1);
  z2 = calculCoordButterFly(3*(A)+2, 3*(B)+2, 3*C+2, 3*D+2, 3*E+2, 3*F+2, 3*G+2, 3*H+2);

  TabVertices[3*TabIndices[3*(oldNbPolys+i)+1]] =   x2;
  TabVertices[3*TabIndices[3*(oldNbPolys+i)+1]+1] = y2;
  TabVertices[3*TabIndices[3*(oldNbPolys+i)+1]+2] = z2;

  TabIndices[3*k]   = TabIndices[3*(oldNbPolys+i)+1]; // Ajout triangle (A,B,MILIEU)
  TabIndices[3*k+1] = A;
  TabIndices[3*k+2] = B;
  k++;

  // troisieme sommet

    A = TabIndices[3*i+2];
    B = TabIndices[3*i+1];
    C = TabIndices[3*i];
    std::cout<<"Point A :"<<A<<" point B : "<<B<<" point C : "<<C<<"\n";

    D = -1;
    G = -1;
    H = -1;
    F = -1;
    E = -1;

    D = indicePointRestantTriangle(A,B, oldNbPolys, C); // Point D, symétrique de C par rapport à AB
    G = indicePointRestantTriangle(D,B, oldNbPolys, A); // Point G, symétrique de A par rapport à DB
    H = indicePointRestantTriangle(A,D, oldNbPolys, B); // Point H, symétrique de B par rapport à AD
    F = indicePointRestantTriangle(A,C, oldNbPolys, B); // Point F, symétrique de B par rapport à AC
    E = indicePointRestantTriangle(C,B, oldNbPolys, A); // Point E, symétrique de A par rapport à CB



  x3 = calculCoordButterFly(3*(A) , 3*(B), 3*C, 3*D, 3*E, 3*F, 3*G, 3*H);
  y3 = calculCoordButterFly(3*(A)+1, 3*(B)+1, 3*C+1, 3*D+1, 3*E+1, 3*F+1, 3*G+1, 3*H+1);
  z3 = calculCoordButterFly(3*(A)+2, 3*(B)+2, 3*C+2, 3*D+2, 3*E+2, 3*F+2, 3*G+2, 3*H+2);

  TabVertices[3*TabIndices[3*(oldNbPolys+i)+2]] =   x3;
  TabVertices[3*TabIndices[3*(oldNbPolys+i)+2]+1] = y3;
  TabVertices[3*TabIndices[3*(oldNbPolys+i)+2]+2] = z3;

  TabIndices[3*k]   = TabIndices[3*(oldNbPolys+i)+2]; // Ajout triangle (A,B,MILIEU)
  TabIndices[3*k+1] = A;
  TabIndices[3*k+2] = B;
  k++;



  // premier triangle
    TabIndices[3*k]   = TabIndices[3*(oldNbPolys+i)];
    TabIndices[3*k+1] = TabIndices[3*(oldNbPolys+i)+1];
    TabIndices[3*k+2] = TabIndices[3*i+1];
/*    std::cout<<"premier triangle"<<std::endl;
    std::cout<< TabIndices[3*k] <<std::endl;
    std::cout<< TabIndices[3*k+1]<<std::endl;
    std::cout<< TabIndices[3*k+2]<<std::endl<<std::endl;*/
    k++;

  // deuxieme triangle

    TabIndices[3*k]   = TabIndices[3*(oldNbPolys+i)];
    TabIndices[3*k+1] = TabIndices[3*i];
    TabIndices[3*k+2] = TabIndices[3*(oldNbPolys+i)+2];
/*    std::cout<<"deuxieme triangle"<<std::endl;
    std::cout<< TabIndices[3*k] <<std::endl;
    std::cout<< TabIndices[3*k+1]<<std::endl;
    std::cout<< TabIndices[3*k+2]<<std::endl<<std::endl;*/
    k++;

  // troisieme triangle
    TabIndices[3*k]   = TabIndices[3*i+2];
    TabIndices[3*k+1] = TabIndices[3*(oldNbPolys+i)+1];
    TabIndices[3*k+2] = TabIndices[3*(oldNbPolys+i)+2];
/*    std::cout<<"troisieme triangle"<<std::endl;
    std::cout<< TabIndices[3*k] <<std::endl;
    std::cout<< TabIndices[3*k+1]<<std::endl;
    std::cout<< TabIndices[3*k+2]<<std::endl<<std::endl;*/
    k++;

}
  // Reajustement des normales
  calculNormalesAuxTriangles(nbPolys);
  // Recalcul des voisins
  fillVoisinTab(TabIndices, TabVoisins, nbPolys);

}