// 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::makerevolution(nbspline & myBSpline)
{
  // you have to implement this
}

void nbsplinesurface::makecoons(const nbspline &Cu0,const nbspline &Cu1,const nbspline &Cv0,const nbspline &Cv1,data_container& data)
{
  nbsplinesurface S1,S2,S3;
  // you have to implement this
}