#ifndef TRIANGLEUNIFORMSUBDIVISION_H
#define TRIANGLEUNIFORMSUBDIVISION_H

#include <vector>
using namespace std;

#include <print.h>
#include <point.h>

/*!
\brief Uniform triangle subdivision on indexed triangles.

The references get where to write the triangles and points too
 when the division functions are called. The are the output.
 The division functions are the input.

This class is non-destructive in the sense that geometry
  is only added to vi and not deleted from it.  So if vi
  already contains geometry this is not destroyed.
*/
template <typename T>
class triangleuniformsubdivision
{
public:

  /** Newly created triangles are appended to vi.*/
  vector< point3<uint> > & vi;
  /** Newly created points are appended to pts. */
  vector< T > & pts;

  /** Pass where to write the new triangles and points too. */
  triangleuniformsubdivision
  ( 
    vector< point3<uint> > & vi_,
    vector< T > & pts_
  ) : vi(vi_), pts(pts_) {}

  /** Divide a triangle one depth into 4 triangles adding 
      the new geometry to vi. */
  void divide( point3<uint> const & ti );

  /** Divide each triangle in v2 one depth into 4 triangles adding 
      the new geometry to vi. */
  void divide( vector< point3<uint> > const & v2 )
    { uintc sz(v2.size()); for (uint i=0; i<sz; ++i) divide(v2[i]); } 

  /** Divide a triangle depth times by dividing one triangle into four adding
      the new geometry to vi. */
  void divide( point3<uint> const & ti, uintc depth );

  /** Divide each triangle in v2 depth times dividing on triangle into 4 triangles adding 
      the new geometry to vi. */
  void divide( vector< point3<uint> > const & v2, uintc depth )
    { uintc sz(v2.size()); for (uint i=0; i<sz; ++i) divide(v2[i],depth); } 

};


//-------------------------------------------------------------
//  Implementation.

template <typename T>
void triangleuniformsubdivision<T>::divide
( 
  point3<uint> const & ti, 
  uintc depth 
)
{
  static vector< point3<uint> > v0;
  static vector< point3<uint> > v1;
  static vector< point3<uint> >* pv[2];

  v0.clear();
  v1.clear();

  pv[0] = & v0;
  pv[1] = & v1;

  static triangleuniformsubdivision<T> v0d(v0,pts);
  static triangleuniformsubdivision<T> v1d(v1,pts);

  static triangleuniformsubdivision<T>* vd[2];
  vd[0] = & v0d;
  vd[1] = & v1d;

  uint current=0;
  pv[1]->push_back(ti);

  uint c2;
  for (uint i=0; i<depth; ++i)
  {
    c2 = (current+1)%2;
    pv[current]->clear();
    vd[current]->divide(*pv[c2]);
    current = (current+1)%2;
  }

  // Point to where the current geometry is.
  current = (current+1)%2;
  uintc sz = pv[current]->size();
  for (uint i=0; i<sz; ++i)
    vi.push_back( (*pv[current])[i] );
}

template <typename T>
void triangleuniformsubdivision<T>::divide( point3<uint> const & ti )
{
  uint i[6];
  i[0] = ti.x;
  i[1] = ti.y;
  i[2] = ti.z;
  i[3] = pts.size(); 
  i[4] = i[3]+1;
  i[5] = i[4]+1;
    
  vi.push_back( point3<uint>(i[3],i[1],i[4]) );
  vi.push_back( point3<uint>(i[5],i[3],i[4]) );
  vi.push_back( point3<uint>(i[5],i[4],i[2]) );
  vi.push_back( point3<uint>(i[0],i[3],i[5]) );
    
  T const & p0( pts[i[0]] );
  T const & p1( pts[i[1]] );
  T const & p2( pts[i[2]] );

// Problem with temporaries, explicit code because
//   compiler generated crap.

  T q0(p0);
  q0 += p1;
  q0 *= 0.5;
  T q1(p1);
  q1 += p2;
  q1 *= 0.5;
  T q2(p2);
  q2 += p0;
  q2 *= 0.5;

  pts.push_back(q0);
  pts.push_back(q1);
  pts.push_back(q2);

/*
cout << SHOW(q0) << endl;
cout << SHOW(q1) << endl;
cout << SHOW(q2) << endl;
*/

    
/*
cout << SHOW(p0) << endl;
cout << SHOW(p1) << endl;
cout << SHOW(p2) << endl;

    pts.push_back((p0+p1)*0.5);
    pts.push_back((p1+p2)*0.5);
    pts.push_back((p2+p0)*0.5);

cout << SHOW(p0+p1) << endl;
cout << SHOW(p1+p2) << endl;
cout << SHOW(p2+p0) << endl;

cout << SHOW(p0) << endl;
cout << SHOW(p1) << endl;
cout << SHOW(p2) << endl;
*/

}








#endif