// Gnurbs - A curve and surface library
// Copyright (C) 2008-2026 Eric Bechet
//
// See the LICENSE file for contributions and license information.
// Please report all bugs and problems to <bechet@cadxfem.org>.
//

#include <iostream>
#include <vector>
#include <deque>

#include "ndisplay.h"
#include "nbsplinesurface.h"
#include "nbeziercurve.h"

void initsurface(nbsplinesurface &surf,double DX,double DY,double DZ); //fonction d'initalisation de la surface
double cos_angle(double V1[2],double M1[2][2],double V2[2]);

struct info_s
{
  npoint3 pt;
  npoint3 dir1;
  npoint3 dir2;
  int info;
};


int main(void)
{
  data_container data;
  ndisplay display;
  double delta=0.005;
  int step=50;
  int limit=500;
  std::vector<npoint3> liste; // la liste des points de la geodesique
  std::vector<npoint3> liste2; // la liste des points de la geodesique
  nbsplinesurface surface(11,3,11,3); // la surface
  initsurface(surface,0,0,0);
  npoint3 point_depart(0.1,0.1,0);

  
  std::deque<std::pair<int,int> > cells;
  const int MAXU=70;
  const int MAXV=70;
  const double cellsize = 0.025;
  const int steps = 3;
  info_s *tab[MAXU][MAXV];
  for (int u=0;u<MAXU;++u)
    for (int v=0;v<MAXV;++v)
      tab[u][v]=0;
  npoint3 dir1(1,0,0);
  npoint3 dir2(1,0,0);
  npoint3 newp1,newdir1,newdirr1;
  npoint3 newp2,newdir2,newdirr2;
  npoint3 newp01,newdir01;
  npoint3 newp02,newdir02;
  
  npoint p;
  double angle=3.1415926536/2;
  double angle1=angle;
  double angle2=angle;

  info_s* nc= new info_s;
  nc->pt=point_depart;
  nc->dir1=dir1;
  surface.rotate(point_depart[0],point_depart[1],dir1,angle,dir2);
  nc->dir2=dir2;
  nc->info=3;
  tab[0][0]=nc;
  newp1=point_depart;
  newdir1=dir1;
  for (int u=1;u<MAXU;++u)
  {
    surface.walk(newp1[0],newp1[1],newdir1,cellsize,newp1,newdir1,steps);
//    newp1.print(cout);
    nc= new info_s;
    nc->pt=newp1;
    nc->dir1=newdir1;
    nc->info=1;
    tab[u][0]=nc;
  }
  newp2=point_depart;
  newdir2=dir2;
  
  for (int v=1;v<MAXV;++v)
  {
    surface.walk(newp2[0],newp2[1],newdir2,cellsize,newp2,newdir2,steps);
    nc= new info_s;
    nc->pt=newp2;
    nc->dir2=newdir2;
    nc->info=2;
    tab[0][v]=nc;
  }
  cells.push_back(std::pair<int,int>(1,1));

  while (!cells.empty())
  {
    int u=cells[0].first;
    int v=cells[0].second;
    cells.pop_front();
    if (v>=MAXV||u>=MAXU) continue;
    int info;
    point_depart=tab[u-1][v-1]->pt;
    newp1=tab[u][v-1]->pt;
    newp2=tab[u-1][v]->pt;
    newdir1=tab[u][v-1]->dir1;
    newdir2=tab[u-1][v]->dir2;
    info=tab[u-1][v-1]->info;
    for(int i=0;i<100;++i)
    {

      if (! surface.walk(newp1[0],newp1[1],newdir2,cellsize,newp01,newdir02,steps)) break;

      if (! surface.walk(newp2[0],newp2[1],newdir1,cellsize,newp02,newdir01,steps)) break;
  
      double M[2][2];

      npoint3 newdirr10b=newp02-newp1;
      npoint3 newdirr10a=newp01-newp1;
      surface.fundamental_forms(newp1[0],newp1[1],M);
      int sign1 = ((newdirr10a[0]*newdirr10b[1]-newdirr10a[1]*newdirr10b[0])>=0)?1:-1;
      double dangle1=acos(cos_angle(newdirr10a.array(),M,newdirr10b.array()))*sign1;

      
      
      npoint3 newdirr20b=newp01-newp2;
      npoint3 newdirr20a=newp02-newp2;
      surface.fundamental_forms(newp2[0],newp2[1],M);
      int sign2 = ((newdirr20a[0]*newdirr20b[1]-newdirr20a[1]*newdirr20b[0])>=0)?1:-1;
      double dangle2=acos(cos_angle(newdirr20a.array(),M,newdirr20b.array()))*sign2;


      if (fabs(dangle1)+fabs(dangle2)<=2e-6) break;
      surface.rotate(newp2[0],newp2[1],newdir1,dangle2,newdir1);
      surface.rotate(newp1[0],newp1[1],newdir2,dangle1,newdir2);      
    }
    nc= new info_s;
    nc->pt=newp01;
    nc->dir1=newdir01;
    nc->dir2=newdir02;
    nc->info=0;
    tab[u][v]=nc;
    tab[u-1][v]->info+=1;
    if (tab[u-1][v]->info==3) cells.push_back(std::pair<int,int>(u,v+1));
    tab[u-1][v]->dir1=newdir1;
    tab[u][v-1]->info+=2;
    if (tab[u][v-1]->info==3) cells.push_back(std::pair<int,int>(u+1,v));
    tab[u][v-1]->dir2=newdir2;
  }

  properties pr=data.getproppoints();
  pr.pointsize=3.0;
  pr.c=color(255,0,0);
  data.setproppoints(pr);
  for (int u=0;u<MAXU-1;++u)
    for (int v=0;v<MAXV-1;++v)
    {
      double M[2][2];
      npoint p,pp,ppp,pppp;
      npoint3 p0,v1,v2;
      npoint3 p1,p2,p3;
      bool ok=true;
      if (tab[u][v])
      {
        p0=tab[u][v]->pt;
        v1=tab[u][v]->dir1;
        v2=tab[u][v]->dir2;
        surface.P(p0[0],p0[1],p);
        surface.fundamental_forms(p0[0],p0[1],M);
        double cangle = cos_angle(v1.array(),M,v2.array());
        if (cangle>=0)
        {
          data.setcolorpoints(color(255*cangle,0,255-255*cangle));
          data.setcolorquads(color(255*cangle,0,255-255*cangle));
        }
        if (cangle<0)
        {
          data.setcolorpoints(color(0,-255*cangle,255+255*cangle));
          data.setcolorquads(color(0,-255*cangle,255+255*cangle));
        }
//        data.add_point(p);
      } else ok=false;
      if (ok)
      {
        if(tab[u+1][v])
          if (tab[u][v+1])
            if (tab[u+1][v+1])
            {
              p1=tab[u+1][v]->pt;
              p2=tab[u+1][v+1]->pt;
              p3=tab[u][v+1]->pt;
              surface.P(p1[0],p1[1],pp);
              surface.P(p2[0],p2[1],ppp);
              surface.P(p3[0],p3[1],pppp);
              data.add_quad(quad(p,pp,ppp,pppp));
            }
      }
    }
  

  data.setcolorquads(color(127,127,127,192));
//  surface.Display(data);
  display.init_data(data);
  display.display();
  return 0;
}


//fonction d'initalisation d'une surface
void initsurface(nbsplinesurface &surf,double DX,double DY,double DZ)
{
  for (int i=0;i<surf.nb_CP_u();++i)
  {
    for (int j=0;j<surf.nb_CP_v();++j)
    {
      double x=i*2/((double)surf.nb_CP_u()-1)-1;
      double y=j*2/((double)surf.nb_CP_v()-1)-1;
      double z=0.;
      if ((i-1)/3==(surf.nb_CP_u()/6)&&((j-1)/3==(surf.nb_CP_v()/6)))
      {
       if ((i==surf.nb_CP_u()/2)&&(j==surf.nb_CP_v()/2)) z=0.0;
           else z=0.4;
      }
      else
        z=0.0;

      npoint p(x+DX,y+DY,z+DZ,1.0);
      surf.CP(i,j)=p;
    }
  }
  int endsu=surf.degree_u()+1;
  int du=surf.degree_u();
  int cpsu=surf.nb_CP_u()-1;
  for (int i=0;i<endsu;++i)
  {
    surf.u(i)=0.0;
    surf.u(cpsu+du+1-i)=(cpsu+du+2-2*endsu+1);
  }
  for (int i=1;i<=(cpsu+du+2-2*endsu);++i)
  {
    surf.u(i+endsu-1)=i;
  }
  int endsv=surf.degree_v()+1;
  int dv=surf.degree_v();
  int cpsv=surf.nb_CP_v()-1;

  for (int i=0;i<endsv;++i)
  {
    surf.v(i)=0.0;
    surf.v(cpsv+dv+1-i)=(cpsv+dv+2-2*endsv+1);
  }
  for (int i=1;i<=(cpsv+dv+2-2*endsv);++i)
  {
    surf.v(i+endsv-1)=i;
  }
}



double cos_angle(double V1[2],double M1[2][2],double V2[2])
{
  double a=0;
  double b=0;
  double c=0;
  for (int i=0;i<2;++i)
  {
    for (int j=0;j<2;++j)
    {
      a+=V1[i]*M1[i][j]*V2[j];
      b+=V1[i]*M1[i][j]*V1[j];
      c+=V2[i]*M1[i][j]*V2[j];
    }
  }
  double sqbc=sqrt(fabs(b*c));
  if (sqbc==0) return 0.;
  double r=a/sqbc;
  if (r>1.0) return 1.0;
  if (r<-1.0) return -1.0;
  return r;
}