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

#include <cmath>
#include <iostream>


int nbsplinesurface::degree(void) const
{
  if (degu>degv) return degu;
  else return degv;
}

int nbsplinesurface::findspan(const std::vector<double> &knots,int nknots,int deg,double u) const
{
  int inter;
  int nCP=nknots-deg-1;
  if (u>=knots[nCP]) return nCP-1;
  if (u<=knots[deg]) return deg;
  //else
  //{
  int low=deg,high=nCP+1;
  int mid= (low+high) /2;
  while ((u<knots[mid]) || (u>= knots[mid+1]))
    //while ((u<knots[mid-1]) || (u>= knots[mid]))
  {
    if (u<knots[mid]) high=mid;
    else low=mid;
    //if (u<knots[mid-1]) high=mid; else low=mid;
    mid= (low+high) /2;
  }
  //inter = mid - 1;
  //}
  return mid;
  //return inter;
}


void nbsplinesurface::findspanmult(const std::vector<double> &knots,int nknots, int deg,double u,int &k,int &s)
{
  int nCP=nknots-deg-1;
  s=0;
  if (u>=knots[nCP])
  {
    k=nCP;
    while (knots[nCP]==knots[k+s])
    {
      s++;
      if (s>deg) break;
    }
    return;
  }
  if (u<=knots[deg])
  {
    k=deg;
    while (knots[deg]==knots[k-s])
    {
      s++;
      if (s>deg) break;
    }
    return;
  }
  int low=deg,high=nCP;
  int mid= (low+high) /2;
  while ((u<knots[mid]) || (u>= knots[mid+1]))
  {
    if (u<knots[mid]) high=mid;
    else low=mid;
    mid= (low+high) /2;
  }
  k=mid;
  while (u==knots[k-s])
  {
    s++;
    if (s>deg) break;
  }
}





void nbsplinesurface::basisfuns(const std::vector<double> &knots,int nknots, int deg, int i, double para, std::vector<double> &N) const
{
  N[0]=1.0;
  std::vector<double> left(deg+1);
  std::vector<double> right(deg+1);
  for (int j=1;j<=deg;++j)
  {
    left[j]=para-knots[i+1-j];
    right[j]=knots[i+j]-para;
    double saved=0.0;
    for (int r=0;r<j;++r)
    {
      double temp=N[r]/ (right[r+1]+left[j-r]);
      N[r]=saved+right[r+1]*temp;
      saved=left[j-r]*temp;
    }
    N[j]=saved;
  }
}




void nbsplinesurface::gradbasisfuns(const std::vector<double> &knots,int nknots,int deg,int i, double u,int n,std::vector<std::vector<double> > &N) const
{
  std::vector<double> left(deg+1);
  std::vector<double> right(deg+1);
  std::vector<std::vector<double> > ndu(deg+1);
  for (int j=0;j<=deg;++j) ndu[j].resize(deg+1);
  std::vector<double> a[2];
  for (int j=0;j<2;++j) a[j].resize(deg+1);
  ndu[0][0]=1.0;
  for (int j=1;j<=deg;++j)
  {
    left[j]=u-knots[i+1-j];
    right[j]=knots[i+j]-u;
    double saved=0.0;
    for (int r=0;r<j;++r)
    {
      ndu[j][r] = right[r+1]+left[j-r];
      double temp=ndu[r][j-1]/ ndu[j][r];
      ndu[r][j]=saved+right[r+1]*temp;
      saved=left[j-r]*temp;
    }
    ndu[j][j]=saved;
  }
  for (int j=0;j<=deg;++j)
  {
    N[0][j]=ndu[j][deg];
  }
  int r;
  for (r=0;r<=deg;++r)
  {
    int s1=0,s2=1;
    a[0][0]=1.0;
    for (int k=1;k<=n;++k)
    {
      double d=0.0;
      int rk=r-k,pk=deg-k;
      if (r>=k)
      {
        a[s2][0]=a[s1][0]/ndu[pk+1][rk];
        d=a[s2][0]*ndu[rk][pk];
      }
      int j,j1,j2;
      if (rk>=-1) j1=1;
      else j1=-rk;
      if (r-1<=pk) j2=k-1;
      else j2=deg-r;
      for (j=j1;j<=j2;++j)
      {
        a[s2][j]=(a[s1][j]-a[s1][j-1])/ndu[pk+1][rk+j];
        d+=a[s2][j]*ndu[rk+j][pk];
      }
      if (r<=pk)
      {
        a[s2][k] = -a[s1][k-1]/ndu[pk+1][r];
        d+=a[s2][k]*ndu[r][pk];
      }
      N[k][r]=d;
      j=s1;s1=s2;s2=j;
    }
  }
  r=deg;
  for (int k=1;k<=n;++k)
  {
    for (int j=0;j<=deg;++j) N[k][j]*=r;
    r*=(deg-k);
  }
}



void nbsplinesurface::curvepoint(double u, double v, npoint &C) const
{
  std::vector<double> Nu(degu+1);
  std::vector<double> Nv(degv+1);
  int spanu=findspan(knotsu,nknotsu,degu,u);
  int spanv=findspan(knotsv,nknotsv,degv,v);
  basisfuns(knotsu,nknotsu,degu,spanu,u,Nu);
  basisfuns(knotsv,nknotsv,degv,spanv,v,Nv);
  for (int i=0;i<4;++i) C[i]=0.;
  for (int i=0;i<=degu;++i)
  {
    int debu=spanu-degu+i;
    for (int j=0;j<=degv;++j)
    {
      int debv=spanv-degv+j;
      C=C+val[debu+nCPu*debv]*Nu[i]*Nv[j];
    }
  }
}


void nbsplinesurface::dercurvepoint(double u, double v, int nu,int nv,npoint &C) const
{
  std::vector<std::vector<double> > Nu(nu+1);
  std::vector<std::vector<double> > Nv(nv+1);
  for (int i=0;i<nu+1;++i) Nu[i].resize(degu+1);
  for (int i=0;i<nv+1;++i) Nv[i].resize(degv+1);
  int spanu=findspan(knotsu,nknotsu,degu,u);
  int spanv=findspan(knotsv,nknotsv,degv,v);
  gradbasisfuns(knotsu,nknotsu,degu,spanu,u,nu,Nu);
  gradbasisfuns(knotsv,nknotsv,degv,spanv,v,nv,Nv);
  for (int i=0;i<4;++i) C[i]=0.;
  for (int i=0;i<=degu;++i)
  {
    int debu=spanu-degu+i;
    for (int j=0;j<=degv;++j)
    {
      int debv=spanv-degv+j;
      C=C+val[debu+nCPu*debv]*Nu[nu][i]*Nv[nv][j];
    }
  }
}

void nbsplinesurface::P(double u,double v, npoint& ret) const
{
  curvepoint(u, v, ret);
}

void nbsplinesurface::dPdu(double u,double v, npoint& ret) const
{
  dercurvepoint(u, v, 1,0,ret);
}

void nbsplinesurface::dPdv(double u,double v, npoint& ret) const
{
  dercurvepoint(u, v,0,1, ret);
}

void nbsplinesurface::d2Pdu2(double u,double v, npoint& ret) const
{
  dercurvepoint(u, v,2,0, ret);
}

void nbsplinesurface::d2Pdv2(double u,double v, npoint& ret) const
{
  dercurvepoint(u, v, 0,2,ret);
}

void nbsplinesurface::d2Pdudv(double u,double v, npoint& ret) const
{
  dercurvepoint(u, v,1,1, ret);
}


void nbsplinesurface::saturate_u()
{
  int curn=degu;
  double curk;
  while (curn< (nknotsu-degu))
  {
    curk=knotsu[curn];
    insertknot_u(curk,degu);
    while ((knotsu[curn]==curk) && (curn< (nknotsu-degu))) curn++;
  }
}

void nbsplinesurface::saturate_v()
{
  int curn=degv;
  double curk;
  while (curn< (nknotsv-degv))
  {
    curk=knotsv[curn];
    insertknot_v(curk,degv);
    while ((knotsv[curn]==curk) && (curn< (nknotsv-degv))) curn++;
  }
}

void nbsplinesurface::insertknot_u(double u,int r)
{
  int k,s;
  findspanmult(knotsu,nknotsu,degu,u,k,s);
  if ((r+s) >degu) r=degu-s;
  if (r>0)
  {
    insertknotu_(u,k,s,r);
  }
}

void nbsplinesurface::insertknot_v(double v,int r)
{
  int k,s;
  findspanmult(knotsv,nknotsv,degv,v,k,s);
  if ((r+s) >degv) r=degv-s;
  if (r>0)
  {
    insertknotv_(v,k,s,r);
  }
}

void nbsplinesurface::insertknotu_(double u,int k,int s,int r)
{
  std::vector<double> newknots(nknotsu+r);
  nbsplinesurface newval(nCPu+r,1,nCPv,1);
  int newnCP=nCPu+r;
  int newnknots=nknotsu+r;
  int mp=nCPu+degu;

  for (int i=0;i<=k;++i) newknots[i]=knotsu[i];
  for (int i=1;i<=r;++i) newknots[i+k]=u;
  for (int i=k+1;i<=mp;++i) newknots[i+r]=knotsu[i];

  for (int v=0;v<nCPv;++v)
  {
    std::vector<npoint> R(nCPu+r);
    for (int i=0;i<=k-degu;++i) newval.CP(i,v) =CP(i,v);
    for (int i=k-s;i<nCPu;++i) newval.CP(i+r,v) =CP(i,v);
    for (int i=0;i<=degu-s;++i) R[i]=CP(k-degu+i,v);
    int L=0;
    for (int j=1;j<=r;++j)
    {
      L=k-degu+j;
      for (int i=0;i<=degu-j-s;++i)
      {
        double alpha= (u-knotsu[L+i]) / (knotsu[i+k+1]-knotsu[L+i]);
        R[i]=R[i+1]*alpha+R[i]* (1.0-alpha);
      }
      newval.CP(L,v) =R[0];
      newval.CP(k+r-j-s,v) =R[degu-j-s];
    }
    for (int i=L+1;i<k-s;++i)
    {
      newval.CP(i,v) =R[i-L];
    }
  }
  knotsu=newknots;
  val=newval.val;
  nCPu=newnCP;
  nknotsu=newnknots;
}


void nbsplinesurface::insertknotv_(double v,int k,int s,int r)
{
  std::vector<double> newknots(nknotsv+r);
  nbsplinesurface newval(nCPu,1,nCPv+r,1);
  int newnCP=nCPv+r;
  int newnknots=nknotsv+r;
  int mp=nCPv+degv;

  for (int i=0;i<=k;++i) newknots[i]=knotsv[i];
  for (int i=1;i<=r;++i) newknots[i+k]=v;
  for (int i=k+1;i<=mp;++i) newknots[i+r]=knotsv[i];

  for (int u=0;u<nCPu;++u)
  {
    std::vector<npoint> R(nCPv+r);
    for (int i=0;i<=k-degv;++i) newval.CP(u,i) =CP(u,i);
    for (int i=k-s;i<nCPv;++i) newval.CP(u,i+r) =CP(u,i);
    for (int i=0;i<=degv-s;++i) R[i]=CP(u,k-degv+i);
    int L=0;
    for (int j=1;j<=r;++j)
    {
      L=k-degv+j;
      for (int i=0;i<=degv-j-s;++i)
      {
        double alpha= (v-knotsv[L+i]) / (knotsv[i+k+1]-knotsv[L+i]);
        R[i]=R[i+1]*alpha+R[i]* (1.0-alpha);
      }
      newval.CP(u,L) =R[0];
      newval.CP(u,k+r-j-s) =R[degv-j-s];
    }
    for (int i=L+1;i<k-s;++i)
    {
      newval.CP(u,i) =R[i-L];
    }
  }
  knotsv=newknots;
  val=newval.val;
  nCPv=newnCP;
  nknotsv=newnknots;
}

void nbsplinesurface::degree_elevation_u()
{
  saturate_u();
  nbsplinesurface buf;
  buf.nCPu= (nb_CP_u() + (nb_CP_u()-1) /degree_u());
  buf.nCPv=nb_CP_v();
  buf.degu=degu+1;
  buf.degv=degv;
  buf.knotsu.reserve(nb_knots_u() + (nb_knots_u()-2) /degree_u());
  buf.knotsv=knotsv;
  buf.nknotsv=nknotsv;
  buf.val.resize((nb_CP_u() + (nb_CP_u()-1) /degree_u()) *nb_CP_v());
  buf.knotsu.resize(0);
//  buf.val.resize(0);
  int i=0;
  for (;i<nb_knots_u()-1;++i)
  {
    buf.knotsu.push_back(knotsu[i]);
    if ((knotsu[i]-knotsu[i+1]) !=0) buf.knotsu.push_back(knotsu[i]);
  }
  buf.knotsu.push_back(knotsu[i]);
  buf.knotsu.push_back(knotsu[i]);

  for (int v=0;v<nb_CP_v();++v)
  {
    i=0;
    int k=0;
    buf.CP(k,v) =CP(i,v);
    ++k;
    ++i;
    while (i<nb_CP_u())
    {
      for (int j=1;j<=degu;++j)
      {
        buf.CP(k,v) =CP(i-1,v) + ((double)(degu+1.-j) / (degu+1.)) * (CP(i,v)-CP(i-1,v));
        ++k;
        ++i;
      }
      buf.CP(k,v) =CP(i-1,v);
      ++k;
    }
  }
  buf.nknotsu=buf.knotsu.size();
  *this=buf;
}


void nbsplinesurface::degree_elevation_v()
{
  saturate_v();
  nbsplinesurface buf;
  buf.nCPu=nb_CP_u();
  buf.nCPv= (nb_CP_v() + (nb_CP_v()-1) /degree_v());
  buf.degv=degv+1;
  buf.degu=degu;
  buf.knotsv.reserve(nb_knots_v() + (nb_knots_v()-2) /degree_v());
  buf.knotsu=knotsu;
  buf.nknotsu=nknotsu;
  buf.val.resize((nb_CP_v() + (nb_CP_v()-1) /degree_v()) *nb_CP_u());
  buf.knotsv.resize(0);
//  buf.val.resize(0);
  int i=0;
  for (;i<nb_knots_v()-1;++i)
  {
    buf.knotsv.push_back(knotsv[i]);
    if ((knotsv[i]-knotsv[i+1]) !=0) buf.knotsv.push_back(knotsv[i]);
  }
  buf.knotsv.push_back(knotsv[i]);
  buf.knotsv.push_back(knotsv[i]);

  for (int u=0;u<nb_CP_u();++u)
  {
    i=0;
    int k=0;
    buf.CP(u,k) =CP(u,i);
    ++k;
    ++i;
    while (i<nb_CP_v())
    {
      for (int j=1;j<=degv;++j)
      {
        buf.CP(u,k) =CP(u,i-1) + ((double)(degv+1.-j) / (degv+1.)) * (CP(u,i)-CP(u,i-1));
        ++k;
        ++i;
      }
      buf.CP(u,k) =CP(u,i-1);
      ++k;
    }
  }
  buf.nknotsv=buf.knotsv.size();
  *this=buf;
}



void nbsplinesurface::makecoons(const nbspline &Cu0,const nbspline &Cu1,const nbspline &Cv0,const nbspline &Cv1,data_container* data)
{
  nbsplinesurface S1,S2,S3;
  S1.knotsu=Cu0.knots;
  S1.nknotsu=Cu0.nknots;
  S1.nknotsv=4;
  S1.knotsv.resize(4);
  S1.knotsv[0]=S1.knotsv[1]=Cv0.min_u();
  S1.knotsv[2]=S1.knotsv[3]=Cv0.max_u();
  S1.val.resize(2*Cu0.nCP);
  S1.nCPu=Cu0.nCP;
  S1.nCPv=2;
  S1.degu=Cu0.deg;
  S1.degv=1;
  for (int i =0;i<Cu0.nCP;++i)
  {
    S1.CP(i,0) =Cu0.CP(i);
    S1.CP(i,1) =Cu1.CP(i);
  }

  S2.knotsv=Cv0.knots;
  S2.nknotsv=Cv0.nknots;
  S2.nknotsu=4;
  S2.knotsu.resize(4);
  S2.knotsu[0]=S2.knotsu[1]=Cu0.min_u();
  S2.knotsu[2]=S2.knotsu[3]=Cu0.max_u();
  S2.val.resize(2*Cv0.nCP);
  S2.nCPv=Cv0.nCP;
  S2.nCPu=2;
  S2.degv=Cv0.deg;
  S2.degu=1;
  for (int i =0;i<Cv0.nb_CP();++i)
  {
    S2.CP(0,i) =Cv0.CP(i);
    S2.CP(1,i) =Cv1.CP(i);
  }

  S3.nknotsu=4;
  S3.knotsu.resize(4);
  S3.knotsu[0]=S3.knotsu[1]=Cu0.min_u();
  S3.knotsu[2]=S3.knotsu[3]=Cu0.max_u();
  S3.nknotsv=4;
  S3.knotsv.resize(4);
  S3.knotsv[0]=S3.knotsv[1]=Cv0.min_u();
  S3.knotsv[2]=S3.knotsv[3]=Cv0.max_u();
  S3.val.resize(4);
  S3.nCPu=2;
  S3.nCPv=2;
  S3.degu=1;
  S3.degv=1;
  Cu0.P(Cu0.min_u(),S3.CP(0,0));
  Cu0.P(Cu0.max_u(),S3.CP(1,0));
  Cu1.P(Cu1.min_u(),S3.CP(0,1));
  Cu1.P(Cu1.max_u(),S3.CP(1,1));
  int degu=Cu0.degree();
  int degv=Cv0.degree();
  int nku=Cu0.nb_knots();
  int nkv=Cv0.nb_knots();
  for (int i=1;i<degu;++i)
  {
    S2.degree_elevation_u();
    S3.degree_elevation_u();
  }
  for (int i=1;i<degv;++i)
  {
    S1.degree_elevation_v();
    S3.degree_elevation_v();
  }
  for (int i=degu+1;i<nku-degu;++i)
  {
    S2.insertknot_u(Cu0.u(i),1);
    S3.insertknot_u(Cu0.u(i),1);
  }
  for (int i=degv+1;i<nkv-degv;++i)
  {
    S1.insertknot_v(Cv0.u(i),1);
    S3.insertknot_v(Cv0.u(i),1);
  }

  *this=S1;
  for (int i=0;i<nb_CP();++i) CP(i) += (S2.CP(i)-S3.CP(i));
  if (data)
  {
    for (int i=0;i<nb_CP();++i)
    {
      S1.CP(i)-=npoint(0,-3,0,0);
      S2.CP(i)-=npoint(0,-6,0,0);
      S3.CP(i)-=npoint(0,-9,0,0);
    }
    data->setcolorquads(color(255,0,0));
    S1.Display(*data);
    data->setcolorquads(color(0,255,0));
    S2.Display(*data);
    data->setcolorquads(color(0,0,255));
    S3.Display(*data);
  }
}