// 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 "nbeziercurve.h"

#include <cmath>
//#include <iostream>

double nbeziercurve::Basis(int which,double u_) const
{
  return binom(which,degree()) * (pow(u_,which) *pow(1-u_,degree()-which));
//  return bezierBasisRecurrence(degree(),which,u_);
}

void nbeziercurve::P(double u_, npoint& ret) const
{
  DeCasteljau(u_,ret);
  /*  for (int k =0 ; k<4;++k)   ret[k]=0;
    for (int i = 0; i <= degree(); ++i)
    {
      double weight=Basis(i,u_);
      ret += val[i] * weight;
    }*/
}

void nbeziercurve::dPdu(double u_,npoint& ret) const
{
  std::vector<npoint> vec(2);
  DeCasteljau(u_,1,vec);
  ret=(vec[1]-vec[0])*degree();
}

void nbeziercurve::d2Pdu2(double u_,npoint& ret) const
{
  std::vector<npoint> vec(3);
  DeCasteljau(u_,2,vec);
  int d=degree();
  ret=(vec[2]+vec[0]-2*vec[1])*d*(d-1);
}



int nbeziercurve::binom(int i,int n) const
{
  if (n<0) return 0;
  if (i==0) return 1;
  if (i==n) return 1;
  if (i<0) return 0;
  if (i>n) return 0;
  int b= binom(i-1,n-1) +binom(i,n-1);
  return b;
}


double nbeziercurve::bezierBasisRecurrence(int degree, int which,double u_) const
{
  if (which>degree) return 0.;
  if (which<0) return 0.;
  if (which==degree) return pow(u_,degree);
  if (which==0) return pow(1-u_,degree);
  return (1-u_) *bezierBasisRecurrence(degree-1,which,u_) +
         u_*bezierBasisRecurrence(degree-1,which-1,u_);
}


void nbeziercurve::DeCasteljau(double u_, npoint & ret) const
{
  double u1_ = 1.-u_;
  if (nCP>2)
  {
    std::vector<npoint> vec(nCP-1);
    for (int i=0;i<nCP-1;++i)
    {
      vec[i]=val[i]*u1_+val[i+1]*u_;
    }
    for (int i=nCP-2;i>1;--i)
    {
      for (int j=0;j<i;++j)
      vec[j]=vec[j]*u1_+vec[j+1]*u_;
    }
    ret=vec[0]*u1_+vec[1]*u_;
    return;
  }
  else if (nCP==2)
  {
    ret=val[0]*u1_+val[1]*u_;
  }
  else
    ret=val[0];
}

void nbeziercurve::DeCasteljau(double u_, int stop,std::vector<npoint>& ret) const
{  
  if (stop>=nCP-1)
  {
    for (int i=0;i<nCP;++i)
    {
      ret[i]=val[i];
    }
    return;
  }
  double u1_=1.-u_;
  if (stop==nCP-2)
  {
    for (int i=0;i<nCP-1;++i)
    {
      ret[i]=val[i]*u1_+val[i+1]*u_;
    }
    return;
  }
  std::vector<npoint> vec(nCP-1);
  for (int i=0;i<nCP-1;++i)
  {
    vec[i]=val[i]*u1_+val[i+1]*u_;
  }
  for (int i=nCP-2;i>stop+1;--i)
  {
    for (int j=0;j<i;++j)
    vec[j]=vec[j]*u1_+vec[j+1]*u_;
  }
  for (int i=0;i<stop+1;++i)
  {
    ret[i]=vec[i]*u1_+vec[i+1]*u_;
  }
}


void nbeziercurve::DeCasteljau(double u_, npoint & ret,nbeziercurve &left,nbeziercurve &right) const
{
  left.add_CP(CP(0));
  right.add_CP(CP(nb_CP()-1));
  if (nb_CP() ==1)
    ret= CP(0);
  else
  {
    nbeziercurve tmp(nb_CP()-1);
    for (int i=0;i<tmp.nb_CP();++i)
      tmp.set_CP(i,CP(i) * (1-u_) +CP(i+1) *u_);
    tmp.DeCasteljau(u_,ret,left,right);
  }
}

void nbeziercurve::cut(double u_, npoint & ret,nbeziercurve &left,nbeziercurve &right) const
{
  left.nCP=0;
  left.val.resize(0);
  right.nCP=0;
  right.val.resize(0);
  DeCasteljau(u_,ret,left,right);    // right se retrouve a l'envers...
  for (int i=0;i<right.nb_CP() / 2 ;++i)
  {
    npoint tmp=right.CP(i);
    right.set_CP(i,right.CP(right.nb_CP()-1-i));
    right.set_CP(right.nb_CP()-1-i,tmp);
  }
}


void nbeziercurve::degree_elevation()
{
  int d=degree();
  add_CP(CP(d));
  for (int i=d;i>0;--i)
  {
    set_CP(i,CP(i-1) * (i/ (d+1.)) + CP(i) * ((d+1-i) / (d+1.)));
  }
}

void nbeziercurve::translate(npoint vec)
{
  for (int i=0;i<nCP;++i) CP(i)=CP(i)+vec;
}

void nbeziercurve::Display(data_container& data) const    // version optimisee
{
  properties p,p1;
  p=data.getproppoints();
  p1=p;
  if (p.nb_uid<2) p.nb_uid=2;
  for (int i=0;i< nb_CP();i++)
  {
    p.uids[1]=i;
    data.setproppoints(p);
    data.add_point(CP(i));
  }
  data.setproppoints(p1);
  std::list<nbeziercurve> curves;
  curves.push_back(*this);
  int deep= (int)(log(16.) /log(2.));
  recursivedivide(curves,deep);
  std::list<nbeziercurve>::iterator it=curves.begin();
  while (it!=curves.end())
  {
    nbeziercurve& current=*it;
    npoint dep;
    npoint arr;
    bool first=true;
    for (int i=0;i<current.nb_CP();++i)
    {
      if (first)
      {
        arr=current.CP(i);
        first=false;
      }
      else
      {
        dep=arr;
        arr=current.CP(i);
        line l;
        l.pts[0]=dep;
        l.pts[1]=arr;
        data.add_line(l);
      }
    }
    it++;
  }
}

void nbeziercurve::recursivedivide(std::list<nbeziercurve> &curves,int deep) const
{
  if (deep)
  {
    std::list<nbeziercurve>::iterator it=curves.begin();
    while (it!=curves.end())
    {
      npoint ret;
      std::list<nbeziercurve>::iterator it2=it;
      it++;
      nbeziercurve left(0);
      nbeziercurve right(0);
      it2->cut(0.5, ret,left,right);
      *it2=right;
      curves.insert(it2,left);
    }
    recursivedivide(curves,deep-1);
  }
}



/*
void saveCPrecur(double val[][2], int degree, const char* filename="cprec.txt")
{
  std::ofstream of(filename,std::ios::app);
  for (int i=0;i<=degree;i++)
  {
    of  << val[i][0] << " " << val[i][1]  << std::endl;
  }
  of.close();
}



void recursivedivide( double val[][2], int degree,int deep,int level)
{
  if (deep==level)
  {
    saveCPrecur(val,degree);
    return;
  }

  double P[100][100][2];
  for (int i=0;i<=degree;++i)
  {
    P[i][0][0]=val[i][0];
    P[i][0][1]=val[i][1];
  }

  for (int j=1;j<=degree;++j)
    for (int i=0;i<=(degree-j);++i)
    {
      P[i][j][0]=0.5*P[i+1][j-1][0]+0.5*P[i][j-1][0];
      P[i][j][1]=0.5*P[i+1][j-1][1]+0.5*P[i][j-1][1];
    }

  double valr[100][2];
  double vall[100][2];
  for (int i=0;i<=degree;++i)
  {
    valr[i][0]=P[0][i][0];
    valr[i][1]=P[0][i][1];
  }
  for (int i=0;i<=degree;++i)
  {
    vall[i][0]=P[i][degree-i][0];
    vall[i][1]=P[i][degree-i][1];
  }
  recursivedivide(vall, degree,deep,level+1);
  recursivedivide(valr, degree,deep,level+1);
}


*/