/**********************************************************************************************
 *nd.cpp
 *
 *Implementation file of a nd visualization library.
 *
 *Josh McCoy
 *October 2003
 *
 **********************************************************************************************/
#include "nd.h"

ND nd;
NDCallList NDlist;

/*void translateNd(GLdouble m[6][6], GLdouble d1, GLdouble d2, GLdouble d3, GLdouble d5)
 *
 *Translates along any of the four axes in the transformation matrix given by the
 *first argument.  The four remaining arguments are the distance of the translations
 *according to each axis.  The matrix argument is changed and holds the results of
 *the transformation.
 */
void translateNd(vector< vector<GLdouble> > & m, const vector<GLdouble> & p)
{
 for(size_t i = 0; i < nd.dimension() && (unsigned int)i < p.size(); ++i)
    m[i][nd.dimension()] = p[i];
}

void glTranslateNd(const vector<GLdouble> & p)
{
  vector< vector<GLdouble> > m = nd.modelviewMatrix();
  translateNd(m, p);
  glSetModelviewMatrixNd(m);
}
/*void rotateNd(GLdouble m[6][6], GLdouble theta, GLint axis)
 *
 *Changes a tranformation matrix to reflect a rotation around an axis.  The
 *first argument is the matrix to be changed.  The second argument is the angle
 *of rotation from the axis in degrees.  The last argument is the axis to rotate
 *around in dimension order (0 rotates around the x-axis, 1 rotates around the y-axis,
 *2 rotates around the z-axis, and 3 rotates around the fourth dimensional axis).
 */
void rotateNd(vector< vector<GLdouble> > & m, GLdouble theta, GLint axis)
{
  vector< vector<GLdouble> > r(nd.dimension()+1, vector<GLdouble>(nd.dimension()+1, 0.0));
  setIDmatrix(r);

  r[(axis)%(nd.dimension()+1)][(axis)%(nd.dimension()+1)] = cos(theta * TO_RADIANS);
  r[(axis)%(nd.dimension()+1)][(axis+1)%(nd.dimension()+1)] = -sin(theta * TO_RADIANS);
  r[(axis+1)%(nd.dimension()+1)][(axis)%(nd.dimension()+1)] = sin(theta * TO_RADIANS);
  r[(axis+1)%(nd.dimension()+1)][(axis+1)%(nd.dimension()+1)] = cos(theta * TO_RADIANS);
  
  multMatrixNd(m, r, m);
}

//2plane rotation
void rotateNd(vector< vector<GLdouble> > & m, GLdouble theta, GLint axis1, GLint axis2)
{
  vector< vector<GLdouble> > r(nd.dimension()+1, vector<GLdouble>(nd.dimension()+1, 0.0));
  setIDmatrix(r);

  r[axis1][axis1] = cos(theta * TO_RADIANS);
  r[axis1][axis2] = -sin(theta * TO_RADIANS);
  r[axis2][axis1] = sin(theta * TO_RADIANS);
  r[axis2][axis2] = cos(theta * TO_RADIANS);
  
  multMatrixNd(m, r, m);
}


//rotates internal modelview matrix around an axis
void glRotateNd(GLdouble theta, GLint axis)
{
  vector< vector<GLdouble> > r;
  r = nd.modelviewMatrix();

  rotateNd(r, theta, axis);

  glSetModelviewMatrixNd(r);
  
}

//rotates internal modelview matrix by twoplanes
void glRotateNd(GLdouble theta, GLint axis1, GLint axis2)
{
  vector< vector<GLdouble> > r;
  r = nd.modelviewMatrix();

  rotateNd(r, theta, axis1, axis2);

  glSetModelviewMatrixNd(r);
  
}

/*void multMatrixNd(GLdouble m1[6][6], GLdouble m2[6][6], GLdouble m3[6][6])
 *
 *Multiples the first matrix argument times the second matrix argument and
 *places the result in the thrid matrix agrument.
 */
void multMatrixNd(const vector< vector<GLdouble> > & m1, 
		  const vector< vector<GLdouble> > & m2, 
		  vector< vector<GLdouble> > & m3)
{
  size_t i=0, j=0, k=0;

  vector< vector<GLdouble> > result(nd.dimension()+1, vector<GLdouble>(nd.dimension()+1,0.0)); 
  
  for(i = 0; i < nd.dimension()+1; ++i){
    for(j = 0; j < nd.dimension()+1; ++j){			            
      for(k = 0; k < nd.dimension()+1; ++k){
	result[i][j] += m1[i][k] * m2[k][j];	
      }
    }
  }
    
  m3 = result;
}

/*void perspectiveDivision(GLdouble p[6])
 *
 *Performs perspective division on the homogenous coordinate agrument.
 *The coordinate is changed and stores the result.
 */
void perspectiveDivision(vector<GLdouble> & p, GLint d)
{
  GLint i;
  for(i=0; i <= d; ++i)
    p[i] = p[i] / p[d];
}

/*void matrixVectorProduct(GLdouble m[6][6], GLdouble p[6])
 *
 *Multiples the matrix agrument by the homogeneous coordinate.
 *The result goes in the coordinate specified by the third argument
 */
void matrixCoordinateProduct(const vector< vector<GLdouble> > & m,
			     const vector<GLdouble> & p,
			     vector<GLdouble> & r)
{
  size_t i, j;
  vector<GLdouble> result(nd.dimension()+1, 0.0);

  for(i = 0; i < nd.dimension()+1; ++i)
    for(j = 0; j < nd.dimension()+1; ++j)
      result[i] += m[i][j] *  p[j];

  r = result;
}
/*vector<GLdouble> matrixCoordinateProduct(const vector< vector<GLdouble> > & m,
 *			     const vector<GLdouble> & p)
 *
 *Multiples the matrix agrument by the homogeneous coordinate.
 *The result goes in the coordinate specified by the third argument
 */
vector<GLdouble> matrixCoordinateProduct(const vector< vector<GLdouble> > & m,
			     const vector<GLdouble> & p)
{
  size_t i, j;
  vector<GLdouble> result(nd.dimension()+1, 0.0);

  for(i = 0; i < nd.dimension()+1; ++i)
    for(j = 0; j < nd.dimension()+1; ++j)
      result[i] += m[i][j] *  p[j];

  return result;
}

/*void printCoordinate(GLdouble p[6])
 *
 *Prints the given coordinate to stdout.
 */
void printCoordinate(const vector<GLdouble> & p)
{
  for(unsigned int i = 0; i < p.size(); ++i)
    cout << p[i] << " ";
  cout << endl;
}

/*void printMatrix(GLdouble m[6][6])
 *
 *Prints the given matrix to stdout.
 */
void printMatrix(const vector< vector <GLdouble> > & m)
{
  cout.precision(4);
  
  for(unsigned int i = 0; i < m.size(); ++i){
    for(unsigned int j = 0; j < m.size(); ++j)
      cout <<  m[i][j] << " ";    
    cout << endl;
  }  
}

void glBeginNd(GLint mode)  //push new CallData on the CallList
{
  NDlist.addtoList(mode);
}

void glEndNd()
{

}


void glVertexNdi(const vector<GLint> & v)
{
  NDlist.addVertex(v);
}


void glVertexNdf(const vector<GLdouble> & v)
{
  NDlist.addVertex(v);
}

void setIDmatrix(vector< vector<GLdouble> > & m)
{
  for(unsigned int i = 0; i < m.size(); ++i)
    for(unsigned int j = 0; j < m.size(); ++j)
      if(i == j)
	m[i][j] = 1.0;
      else
	m[i][j] = 0.0;
}


vector< vector<GLdouble> > getIDmatrix(GLint size)
{
  vector< vector<GLdouble> > id(size, vector<GLdouble>(size, 0.0));
  GLint i;
  for(i=0; i < size; ++i)
    id[i][i] = 1.0;

  return id;
}

//returns an identity matrix with size of n+1 x n+1
vector< vector<GLdouble> > getIDmatrix()
{
  return getIDmatrix(nd.dimension()+1);
}

void glDimensions(GLint n)
{
  nd.setDimension(n);
}

//mutator access to internal matrices
void glSetModelviewMatrixNd(const vector< vector<GLdouble> > & m)
{
  nd.setModelviewMatrix(m);
}

void glSetProjectionMatrixNd(const vector< vector<GLdouble> > & m)
{
  nd.setProjectionMatrix(m);
}

void glFrustumNd(const vector<GLdouble> & f)
{
  nd.setFrustum(f);
}

/*************************************************************
 *Implementation of VertexInfo Class
 ************************************************************/

VertexInfo::VertexInfo(vector<GLdouble> flist)
{
  //if the size of the incoming vector isn't one greater than the number
  //of dimensions, pack up to the nth element with 0's and pack the
  //n+1 element with 1.
  unsigned int i;
  _fvert = flist;
  if(_fvert.size() < nd.dimension()+1){
    for(i=_fvert.size(); i < nd.dimension(); ++i)
      _fvert.push_back(2.0);
    _fvert.push_back(1.0);
  }
  _float = true;
  _int = false;
  _color[0]=_color[1]=_color[2]=_color[3]=0.0;
  glGetFloatv(GL_CURRENT_COLOR, _color);
}

VertexInfo::VertexInfo(vector<GLint> ilist)
{
  //if the size of the incoming vector isn't one greater than the number
  //of dimensions, pack up to the nth element with 0's and pack the
  //n+1 element with 1.
  unsigned int i;
  _ivert = ilist;
  if(_ivert.size() < (unsigned int)nd.dimension()+1){
    for(i=_ivert.size(); i < (unsigned int)nd.dimension(); ++i)
      _ivert.push_back(2);
    _ivert.push_back(1);
  }

  _float = false;
  _int = true;
  glGetFloatv(GL_CURRENT_COLOR, _color);
}

void VertexInfo::perspectiveDivisionOnVertex()
{
  perspectiveDivision(_fvert, nd.dimension());
}

void VertexInfo::projectionOnVertex(const vector< vector<GLdouble> > & proj)
{
  matrixCoordinateProduct(proj, _fvert, _fvert);
}


/************************************************************
 *Implementation of CallData Class
 **********************************************************/

CallData::CallData(GLint mode)
{
  _mode = mode;
}

void CallData::addVertex(vector<GLdouble> v)
{
  _vinfo.push_back(VertexInfo(v));
}

void CallData::addVertex(vector<GLint> v)
{
  _vinfo.push_back(VertexInfo(v));
}

void CallData::perspectiveDivisionOnVertices()
{
  unsigned int i;
  for(i = 0; i < _vinfo.size(); ++i)
    _vinfo[i].perspectiveDivisionOnVertex();
}

void CallData::projectionOnVertices(const vector< vector<GLdouble> > & proj)
{
  unsigned int i;
  for(i = 0; i < _vinfo.size(); ++i)
    _vinfo[i].projectionOnVertex(proj);
}

/*************************************************************
 *Implementation of the NDCallList Class
 ***********************************************************/

void NDCallList::addtoList(GLint mode)
{
  _data.push_back(CallData(mode));
}


void NDCallList::toOpenGL3d()
{
  unsigned int modes, vertices;
  vector< vector<GLdouble> > proj;
  vector<GLdouble> frustum = nd.frustum();

  for(; nd.dimension() >= 3; nd.setDimension(nd.dimension()-1)){
  //prepare projection matrix for the new dimension
    proj = nd.projectionMatrix();
    proj[nd.dimension()][nd.dimension()-1]= 1.0 / frustum[2*(nd.dimension()-1)];
    proj[nd.dimension()][nd.dimension()]=0.0;

    printMatrix(proj);
  //translate all coordinates by the projection matrix
    NDlist.projectionOnList(proj);
  //perform perspective division on all coordinates
    NDlist.perspectiveDivisionOnList();
    cout << "----------------" << nd.dimension() << "--------------------" << endl;
    //    printMatrix(nd.projectionMatrix());
    NDlist.printNd();

  }

  for(modes = 0; modes < _data.size(); ++modes){

    glBegin(_data[modes].mode());

    for(vertices = 0; vertices < _data[modes].nVertices(); ++vertices){

      glColor4fv(_data[modes][vertices].color());

      if(_data[modes][vertices].isfloat()){
	vector<GLdouble> v = _data[modes][vertices].flist();
	glVertex3f(v[0], v[1], v[2]);

      }else{
	vector<GLint> v = _data[modes][vertices].ilist();
	glVertex3f(v[0], v[1], v[2]);
      }	
    }

    glEnd();

  }  
}


//need to have vertex output fixed
void NDCallList::printNd()
{
  unsigned int modes, vertices, i;

  for(modes = 0; modes < _data.size(); ++modes){
    
    //glBegin(_data[modes].mode());
    cout << "glBegin(" << _data[modes].mode() << ");" << endl;

    for(vertices = 0; vertices < _data[modes].nVertices(); ++vertices){

      GLfloat *temp;
      temp = _data[modes][vertices].color();
      //      glColor4fv(_data[modes][vertices].color());
      cout << "glColor4fv(" << temp[0] << ", " << temp[1] << ", "
	   << temp[2] << ", " << temp[3] << ");" << endl;
      

      if(_data[modes][vertices].isfloat()){
	vector<GLdouble> v = _data[modes][vertices].flist();
	//glVertex3f(v[0], v[1], v[2]);
	cout << "glVertexNd(";
	for(i=0; i < v.size()-1; ++i)
	  cout << v[i] << ", ";
      cout << v[i] << ");" << endl;
      }else{
	vector<GLint> v = _data[modes][vertices].ilist();
	//glVertex3i(v[0], v[1], v[2]);
	cout << "glVertex3i(" << v[0] << ", " << v[1] << ", "
	     << v[2] << ");" << endl;
      }	
    }
    //glEnd();
    cout << "glEnd();" << endl;

  } 
}

void NDCallList::printOpenGL3d()
{
  unsigned int modes, vertices;

  vector< vector<GLdouble> > proj;
  vector<GLdouble> frustum = nd.frustum();

  for(; nd.dimension() >= 3; nd.setDimension(nd.dimension()-1)){
    //prepare projection matrix for the new dimension
    proj = nd.projectionMatrix();
    proj[nd.dimension()][nd.dimension()-1]= 1.0 / frustum[2*(nd.dimension()-1)];
    proj[nd.dimension()][nd.dimension()]=0.0;
    
    //translate all coordinates by the projection matrix
    NDlist.projectionOnList(proj);
  //perform perspective division on all coordinates
    NDlist.perspectiveDivisionOnList();
    cout << "----------------" << nd.dimension() << "--------------------" << endl;
    printMatrix(proj);
    NDlist.printNd();
    
  }

  cout << "--------------OpenGL Code----------------" << endl;
  
  for(modes = 0; modes < _data.size(); ++modes){

    cout << "glBegin(" << _data[modes].mode() << ");" << endl;

    for(vertices = 0; vertices < _data[modes].nVertices(); ++vertices){

      GLfloat *temp;
      temp = _data[modes][vertices].color();
      cout << "glColor4fv(" << temp[0] << ", " << temp[1] << ", " 
	   << temp[2] << ", " << temp[3] << ");" << endl;
      if(_data[modes][vertices].isfloat()){
	vector<GLdouble> v = _data[modes][vertices].flist();
	cout << "glVertex3f(" << v[0] << ", " << v[1] << ", "
	     << v[2] << ");" << endl;

      }else{
	vector<GLint> v = _data[modes][vertices].ilist();
	cout << "glVertex3i(" << v[0] << ", " << v[1] << ", "
	     << v[2] << ");" << endl;
      }	
    }
    cout << "glEnd();" << endl;
  }  
}


void NDCallList::addVertex(vector<GLdouble> v)
{
  //failure if size of _data is 0
  _data[_data.size()-1].addVertex(v);
}

void NDCallList::addVertex(vector<GLint> v)
{
  //failure if size of _data is 0
  _data[_data.size()-1].addVertex(v);
}

void NDCallList::perspectiveDivisionOnList()
{
  unsigned int i;
  for(i = 0; i < _data.size(); ++i)
    _data[i].perspectiveDivisionOnVertices();
}

void NDCallList::projectionOnList(const vector< vector<GLdouble> > & proj)
{
  unsigned int i;
  for(i = 0; i < _data.size(); ++i)
    _data[i].projectionOnVertices(proj);
}


/*****************************************************************
 *Implementation of the ND Class
 ****************************************************************/

void ND::setDimension(GLint d)
{
  if(d != _d){
    _d = d;
    _modelview.resize(d+1, vector<GLdouble> (d+1, 0.0));
    _projection.resize(d+1, vector<GLdouble> (d+1, 0.0));
    _frustum.resize(d*2, 1.0);
  }
}