#include "loft.h"

namespace other{

double p_len(Array<const TV>& p){
  double l=0;
  int n = (int)p.size();
  for(int i=n-1, j=0; j<n; i=j++){
    l += (p[i]-p[j]).magnitude();
  }
  return l;
}

double p_len_array(Array<const TV>& p, const int start, const int stride){
  double l=0;
  int n = (int)start+stride;
  for(int i=n-1, j=start; j<n; i=j++){
    l += (p[i]-p[j]).magnitude();
  }
  return l;
}


double r_len(Array<const TV>& r){
  double l=0;
  int n = (int)r.size();
  for(int i=0; i<n-1; i++){
    l += (r[i]-r[i+1]).magnitude();
  }
  return l;
}

const Array<Vector<int,3> > loft_tris(Array<const TV> p1, Array<const TV> p2){
  int n = p1.size();
  int m = p2.size();
  double s = p_len(p1)/p_len(p2);
  Array<Vector<int,3> > tris(n+m,false);
  double d=0.;
  int i=0;
  int j=0;
  while(i<n && j<m){
    double v1 = (p1[(i+1)%n]-p1[i]).magnitude();
    double v2 = (p2[(j+1)%m]-p2[j]).magnitude();
    if(v1+d < s*v2){
      tris[i+j].set(i,(i+1)%n,j+n);
      i++;
      d=d+v1;
    }
    else if(v1+d > s*v2){
      tris[i+j].set(i,(j+1)%m + n, j+n);
      j++;
      d=d - s*v2;
    }
    else{
      tris[i+j].set(i, (i+1)%n, (j+1)%m + n);
      tris[i+j+1].set(i, (j+1)%m + n, j+n);
      i++;
      j++;
    }
  }
  tris[i+j].set((i<n)?i:j+n,0,n);
  return tris;
}

//compute the loft mesh, assuming the vertices are {p1,p2}
const Ref<TriangleMesh> compute_loft(Array<const TV> p1, Array<const TV> p2){
  return new_<TriangleMesh>(loft_tris(p1,p2));
}

//compute lofts on an array with stride st, connecting last to first
const Ref<TriangleMesh> compute_lofts_circ(Array<const TV> p, const int st, bool align){
  int n = p.size();
  Array<Vector<int,3> > tris(2*n,false);
  int m = n/st;
  for(int l=0; l<m; ++l){
    Array<TV> p1(st,false);
    Array<TV> p2(st,false);

    int offset = 0;
    if(align){
      real dist = std::numeric_limits<T>::max();
      for(int i=0; i<st; ++i){
        real d = magnitude(p[l*st]-p[i+((l+1)%m)*st]);
        if(d<dist){
          dist = d;
          offset = i;
        }
      }
      std::cout << offset << std::endl;
    }
    for(int i=0; i<st; ++i){
      p1[i] = p[i+l*st];
      p2[i] = p[(i+offset)%st + ((l+1)%m)*st];
    }
    const Array<Vector<int,3> >& temp = loft_tris(p1,p2);
    for(int i=0; i<2*st; ++i){
      int x = (temp[i].x+st*l)%n;
      int y = (temp[i].y+st*l)%n;
      int z = (temp[i].z+st*l)%n;
      tris[i+l*2*st].set(x,y,z);
    }
  }
  return new_<TriangleMesh>(tris);
}



//compute the rail mesh, assuming the vertices are {p1,p2}
const Ref<TriangleMesh> compute_rail(Array<const TV> r1, Array<const TV> r2){
  int n = r1.size();
  int m = r2.size();
  double s = r_len(r1)/r_len(r2);
  Array<Vector<int,3> > tris(n+m-2,false);
  double d=0.;
  int i=0;
  int j=0;
  while(i<n-1 && j<m-1){
    double v1 = (r1[i+1]-r1[i]).magnitude();
    double v2 = (r2[j+1]-r2[j]).magnitude();
    if(v1+d < s*v2){
      tris[i+j].set(i,i+1,j+n);
      i++;
      d=d+v1;
    }
    else if(v1+d > s*v2){
      tris[i+j].set(i,j+1+n, j+n);
      j++;
      d=d - s*v2;
    }
    else{
      tris[i+j].set(i, i+1, j+1+n);
      tris[i+j+1].set(i, j+1+n, j+n);
      i++;
      j++;
    }
  }
  tris[i+j].set((i<n-1)?i:j+n, n-1, n+m-1);
  return new_<TriangleMesh>(tris);
}

}

using namespace other;

void wrap_loft() {
  OTHER_FUNCTION(compute_loft)
  OTHER_FUNCTION(compute_lofts_circ)
  OTHER_FUNCTION(compute_rail)
}