// 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 "nbezierconvexpoly.h"
#include <iostream>

polybary::polybary(int ns_) : vi(ns_)
{
  for (int i=0;i<ns_;++i)
    vi[i]=npoint2(cos(i*2*M_PI/ns_),sin(i*2*M_PI/ns_));
}

void polybary::get_bary(npoint2 v, std::vector<double> &bary) const
{
  int sz=vi.size();
  std::vector<npoint2> si(vi);
  for (int i=0;i<sz;++i)
  {
    si[i]-=v;
    bary[i]=0.0;
  }
  std::vector<double> ri(sz);
  std::vector<double> Ai(sz);
  std::vector<double> Di(sz);
  for (int i=0;i<sz;++i)
  {
    ri[i]=si[i].norm();
    int next=(i+1)%sz;
    Ai[i]=crossprod(si[i],si[next]);
    Di[i]=si[i].dotprod(si[next]);
    if (ri[i]==0.0)
    {
      bary[i]=1.0;
      return;
    }
    if ((Ai[i]==0.0) && (Di[i]<0.0))
    {
      ri[next]=si[next].norm();
      double sum = ri[i]+ri[next];
      bary[i]=ri[next]/sum;
      bary[next]=ri[i]/sum;
      return;
    }
  }

  double W=0.;
  for (int i=0;i<sz;++i)
  {
    int prev;
    if (i)
      prev=i-1;
    else
      prev=sz-1;
    int next=(i+1)%sz;
    if (Ai[prev]!=0.0)
      bary[i]+=(ri[prev]-Di[prev]/ri[i])/Ai[prev];
    if (Ai[i]!=0.0)
      bary[i]+=(ri[next]-Di[i]/ri[i])/Ai[i];
    W+=bary[i];
  }
  for (int i=0;i<sz;++i)
  {
    bary[i]/=W;
  }
}

int binomial_coefficient(int n, int k)
{
    if (0 == k || n == k)
        return 1;
    if (k > n)
        return 0;
    if (k > (n - k))
        k = n - k;
    if (1 == k)
        return n;
    return binomial_coefficient(n-1,k-1)+binomial_coefficient(n-1,k);
}

polyindex::polyindex(const std::vector<int> &vec) : ns(vec.size())
{
  deg=0;
  for (int i=0;i<vec.size();++i) deg+=vec[i];
  tab.resize(deg);
  int k=deg;
  for (int i=0;i<vec.size();++i)
    for (int j=0;j<vec[i];++j)
      tab[--k]=i;
}

int polyindex::get1index(void) const// get an equivalent unique index for storage, alas not dense
{
  if (ns<0) return -1;
  int ret=0;
  for (int i=deg-1;i>=0;--i) ret+=ret*ns+tab[i];
  return ret;
}

void polyindex::print(void) const
{
  std::vector<int> ndx(ns,0);
  for (int i=0;i<deg;++i) { std::cout << tab[i] << " " ; if (tab[i]<ns) ndx[tab[i]]++;}
  std::cout << "|| ";
  for (int i=0;i<ns;++i) { std::cout << ndx[i] << " " ;}
  std::cout << "|| " << get1index();
  std::cout << std::endl;
}

void polyindex::getindex(std::vector< int >& vec) const
{
  if (ns>0)
  {
    vec.resize(ns);std::fill(vec.begin(), vec.end(), 0);
    for (int i=0;i<deg;++i) if ((tab[i]<ns)&&(tab[i]>=0)) vec[tab[i]]++;
  }
}

polyindex& polyindex::operator++() // pre increment
{
  inc();
  return *this;
}

polyindex polyindex::operator++(int) // post increment
{
    polyindex tmp(*this); // copy
    operator++(); // pre-increment
    return tmp;   // return old value
}

bool polyindex::end(void) // tells if the intex went beyond the end.
{
  if (ns<0) return true;
  return false;
}

bool polyindex::inc(void) // increment & tells if it overflows.
{
  int k=0;
  if (deg==0)
  {
    ns=-1;
    return false;
  }
  if (ns<0) return false;
  bool fini=false;
  do
  {
    tab[k]++;
    if (tab[k]>=ns)
    {
      tab[k]=0;
      k++;
      if (k==deg)
      {
        ns=-1;
        return false;
      }
    }
    else
    {
      fini=true;
      for (int kk=0;kk<k;++kk) tab[kk]=tab[k];
    }
  } while ((k<deg)&&(!fini));
  return true;
}

void polyindex::zero(int ns_) // cancel all indexes
{
  for (int i=0;i<deg;++i) tab[i]=0;
  ns=ns_;
}

int polyindex::getnumindexes(void) const
{
  if (ns<0) return 0;
  return binomial_coefficient(deg+ns-1,deg);
}

polyindex polyindex::operator +(const polyindex& other) const// adds two polyindex of same number of ns.
{
  if (other.ns==ns)
  {
    std::vector<int> ndx1(ns),ndx2(ns);
    getindex(ndx1);
    other.getindex(ndx2);
    for (int i=0;i<ns;++i)
      ndx1[i]+=ndx2[i];
    return polyindex(ndx1);
  }
  else
    return polyindex(0,-1);
}

nbezierconvexpoly::nbezierconvexpoly(int degree_, int ns_) : deg(degree_),ns(ns_)
{
  //inspect all points
  polyindex pi(deg,ns);
  nCP=pi.getnumindexes();
  map.reserve(nCP);
  val.reserve(nCP);
  do
  {
    int ndx1=pi.get1index();
    map.push_back(ndx1);
    val.push_back(npoint(0,0,0,1));
  }
  while (pi.inc());
}

void nbezierconvexpoly::P(double u,double v, npoint& ret) const 
{
  polybary pb(ns);
  std::vector<double> coords(ns);
  pb.get_bary(npoint2(u,v),coords);
  getpoint(coords,ret);
}

template<class It, class T >
It binary_find(It first, It last, const T& value)
{
  first = std::lower_bound(first, last, value);
  return ((first != last) && !(value< *first)) ? first : last;
}


const npoint& nbezierconvexpoly::CP(const polyindex &pi) const
{
  return val[std::distance(map.begin(), binary_find(map.begin(), map.end(), pi.get1index()))];
}

npoint nbezierconvexpoly::CP(int which) const
{
  return val[which];
}

npoint& nbezierconvexpoly::CP(int which)
{
    return val[which];
}

npoint& nbezierconvexpoly::CP(const polyindex &pi)
{
  return val[std::distance(map.begin(), binary_find(map.begin(), map.end(), pi.get1index()))];
}

void nbezierconvexpoly::set_CP(const polyindex &pi,const npoint& pt)
{
  val[std::distance(map.begin(), binary_find(map.begin(), map.end(), pi.get1index()))]=pt;
}

void nbezierconvexpoly::set_CP(int which,const npoint& pt)
{
  val[which]=pt;
}

void nbezierconvexpoly::getpoint(std::vector<double> coord, npoint &C) const
{
  if (deg==0) C=val[0];
  else
  {
    nbezierconvexpoly low(deg-1,ns);
    for(polyindex pi(deg-1,ns);!pi.end();++pi)
    {
      npoint pt(0.,0.,0.,0.0);
      int i=0;
      for (polyindex pi1(1,ns);!pi1.end();++pi1)
      {
        pt+=CP(pi+pi1)*coord[i];
        ++i;
      }
      low.CP(pi)=pt;
    }
    low.getpoint(coord,C);
  }
}



void nbezierconvexpoly::Display(data_container& data) const    // default : uses P(u,res)
{
  polybary pb(ns);
  for (int i=0;i< nb_CP();i++)
    data.add_point(CP(i));
  const int nbpts=1*nb_CP() +10;
  for (int i=0;i<pb.vi.size();++i)
  {
    int cur=i;
    int nxt=(i+1)%pb.vi.size();
    npoint2 c(0.0,0.0); // center of the ref element
    npoint2 p1=pb.vi[cur];
    npoint2 p2=pb.vi[nxt];
    std::vector<npoint> dep(nbpts+1);
    std::vector<npoint> arr(nbpts+1);
    bool first=true;
    for (int iu=0;iu<nbpts;++iu)
    {
      npoint2 p3=(p1*iu+p2*(nbpts-iu-1))/(nbpts-1);
      if (first)
      {
        for (int iv=0;iv<nbpts;++iv)
        {
          npoint2 pt=(p3*iv+c*(nbpts-iv-1))/(nbpts-1);
          P(pt[0],pt[1],arr[iv]);
        }
        first=false;
      }
      else
      {
        dep=arr;
        for (int iv=0;iv<nbpts;++iv)
        {
          npoint2 pt=(p3*iv+c*(nbpts-iv-1))/(nbpts-1);
          P(pt[0],pt[1],arr[iv]);
        }
        quad q;
        for (int j=0;j< (nbpts-1);++j)
        {
          q.pts[0]=dep[j];
          q.pts[1]=arr[j];
          q.pts[2]=arr[j+1];
          q.pts[3]=dep[j+1];
          data.add_quad(q);
        }
      }
    }
  }
}