#ifndef PRIMSHPCENTERS_H
#define PRIMSHPCENTERS_H


#include <tetrahedron.h>

#include <gobj.h>
#include <mathlib.h>


//  Display the triangles centers and mesh.
//  Configure by setting the boolean flags.
//  Then call eval() to post the graphics.

#define SQUARED(X) ((X)*(X))



template< typename T >
class visualize_tetrahedron
{
  gobjContainer & gX;
  tetrahedron<T> const & t;
public:

  // Display the triangle.
  bool mesh_on;
  // Display the midpoints.
  bool midpoints_on;
  // Display the centroid. 
  bool midpoint_center_on;
  // Display the altitudes to there opposite points, 
  // and the orthocenter.
  bool altitude_to_point_on;
  // Display the circumcenter.
  bool altitude_on;

  bool innersphere_on;

  bool outersphere_on;

  // Constructor with basic defaults.
  visualize_tetrahedron
  (
    gobjContainer & _gX,
    tetrahedron<T> const & _t
  );

  // Evaluate or send the graphics to
  void eval() const;

};



// --------------------------------------------------------- 
// Implementation
//



template< typename T >
visualize_tetrahedron<T>::visualize_tetrahedron
(
  gobjContainer & _gX,
  tetrahedron<T> const & _t
)
  : gX(_gX), t(_t), mesh_on(true), midpoints_on(true),
    midpoint_center_on(true), altitude_to_point_on(true),
    altitude_on(true), innersphere_on(true), outersphere_on(true)
{
}


template< typename T >
void visualize_tetrahedron<T>::eval() const 
{
  gX.push_back( new gobjglDisable(GL_LIGHTING) );

  if (mesh_on)
  {
  gX.push_back( new gobjglColor3f(0.0,0.0,1.0) );

  gX.push_back( new gobjglBegin(GL_LINES) );

  gX.push_back( new gobjglVertex3f(t.pi[0]) );
  gX.push_back( new gobjglVertex3f(t.pi[1]) );
  gX.push_back( new gobjglVertex3f(t.pi[0]) );
  gX.push_back( new gobjglVertex3f(t.pi[2]) );
  gX.push_back( new gobjglVertex3f(t.pi[0]) );
  gX.push_back( new gobjglVertex3f(t.pi[3]) );
  gX.push_back( new gobjglVertex3f(t.pi[1]) );
  gX.push_back( new gobjglVertex3f(t.pi[2]) );
  gX.push_back( new gobjglVertex3f(t.pi[1]) );
  gX.push_back( new gobjglVertex3f(t.pi[3]) );
  gX.push_back( new gobjglVertex3f(t.pi[2]) );
  gX.push_back( new gobjglVertex3f(t.pi[3]) );

  gX.push_back( new gobjglEnd() );
  }

  point3<double> m0 = t.midpoint(0);
  point3<double> m1 = t.midpoint(1);
  point3<double> m2 = t.midpoint(2);
  point3<double> m3 = t.midpoint(3);


  gobjQuadric* gd = new gobjQuadric();
  gd->radius=.02;
  gX.push_back( gd );

  if (outersphere_on)
  {
  gX.push_back( new gobjglColor3ub(255,20,147) );
  point3<double> c1 = t.circumcenter();
  gX.push_back( new gobjMySphereDraw(c1,gd) );
for (uint i=0; i<4; ++i)
  cout << (c1-t.pi[i]).distance() << endl;
  }

  if (midpoints_on)
  {
  gX.push_back( new gobjglColor3ub(255,0,255) );
  gX.push_back( new gobjMySphereDraw(m0,gd) );
  gX.push_back( new gobjMySphereDraw(m1,gd) );
  gX.push_back( new gobjMySphereDraw(m2,gd) );
  gX.push_back( new gobjMySphereDraw(m3,gd) );
  }

  if (midpoint_center_on)
  {
  gX.push_back( new gobjglColor3ub(0,255,0) );
  gX.push_back( new gobjglBegin(GL_LINES) );
  gX.push_back( new gobjglVertex3f(t.pi[0]) );
  gX.push_back( new gobjglVertex3f(m0) );
  gX.push_back( new gobjglVertex3f(t.pi[1]) );
  gX.push_back( new gobjglVertex3f(m1) );
  gX.push_back( new gobjglVertex3f(t.pi[2]) );
  gX.push_back( new gobjglVertex3f(m2) );
  gX.push_back( new gobjglVertex3f(t.pi[3]) );
  gX.push_back( new gobjglVertex3f(m3) );
  gX.push_back( new gobjglEnd() );
  }


  point3<double> alt0 = t.verticiespoint(0);
  point3<double> alt1 = t.verticiespoint(1);
  point3<double> alt2 = t.verticiespoint(2);
  point3<double> alt3 = t.verticiespoint(3);

  if (altitude_to_point_on)
  {
  gX.push_back( new gobjglColor3ub(255,215,0) );
  //gX.push_back( new gobjMySphereDraw(alt0,gd) );
  gX.push_back( new gobjMySphereDraw(alt1,gd) );
  gX.push_back( new gobjMySphereDraw(alt2,gd) );
  gX.push_back( new gobjMySphereDraw(alt3,gd) );



  gX.push_back( new gobjglBegin(GL_LINES) );
  gX.push_back( new gobjglVertex3f(t.pi[0]) );
  gX.push_back( new gobjglVertex3f(alt0) );
  gX.push_back( new gobjglVertex3f(t.pi[1]) );
  gX.push_back( new gobjglVertex3f(alt1) );
  gX.push_back( new gobjglVertex3f(t.pi[2]) );
  gX.push_back( new gobjglVertex3f(alt2) );
  gX.push_back( new gobjglVertex3f(t.pi[3]) );
  gX.push_back( new gobjglVertex3f(alt3) );
  gX.push_back( new gobjglEnd() );
  }

/*

  if (innersphere_on)
  {
  point3<double> b0 = t.bisectangle(0);
  point3<double> b1 = t.bisectangle(1);
  point3<double> b2 = t.bisectangle(2);
  point3<double> b3 = t.bisectangle(3);
  gX.push_back( new gobjglColor3ub(127,255,212) );
  gX.push_back( new gobjMySphereDraw(b0,gd) );
  gX.push_back( new gobjMySphereDraw(b1,gd) );
  gX.push_back( new gobjMySphereDraw(b2,gd) );
  gX.push_back( new gobjMySphereDraw(b3,gd) );



//  point2<double> x;
//  double r;
//  t.innercircle(x,r);
//  gX.push_back( new gobjMySphereDraw(x,gd) );
//  gX.push_back( new gobjMyCircleDraw(r,point3<double>(x),*gc) );



  gX.push_back( new gobjglBegin(GL_LINES) );
  gX.push_back( new gobjglVertex3f(t.pi[0]) );
  gX.push_back( new gobjglVertex3f(b0) );

  gX.push_back( new gobjglVertex3f(t.pi[1]) );
  gX.push_back( new gobjglVertex3f(b1) );
  gX.push_back( new gobjglVertex3f(t.pi[2]) );
  gX.push_back( new gobjglVertex3f(b2) );
  gX.push_back( new gobjglVertex3f(t.pi[3]) );
  gX.push_back( new gobjglVertex3f(b3) );

  gX.push_back( new gobjglEnd() );
  }
*/

/*
  {
    double a[12];

//    a[0] = t.pi[1].x - t.pi[2].x;  a[1] = t.pi[1].y - t.pi[2].y;  a[2] = t.pi[1].z - t.pi[2].z; 
//    a[3] = (t.pi[1].dot() - t.pi[2].dot())*0.5;
//    a[4] = t.pi[0].x - t.pi[2].x;  a[5] = t.pi[0].y - t.pi[2].y;  a[6] = t.pi[0].z - t.pi[2].z; 
//    a[7] = (t.pi[0].dot() - t.pi[2].dot())*0.5;
//    a[8] = t.pi[0].x - t.pi[1].x;  a[9] = t.pi[0].y - t.pi[1].y;  a[10] = t.pi[0].z - t.pi[1].z; 
//    a[11] = (t.pi[0].dot() - t.pi[1].dot())*0.5;

    a[0] = t.pi[0].x; a[1] = t.pi[0].y; a[2] = t.pi[0].z; a[3] = t.pi[0].dot()*0.5;
    a[4] = t.pi[1].x; a[5] = t.pi[1].y; a[6] = t.pi[1].z; a[7] = t.pi[1].dot()*0.5;
    a[8] = t.pi[2].x; a[9] = t.pi[2].y; a[10] = t.pi[2].z; a[11] = t.pi[2].dot()*0.5;

    gausselim<double,3> g(a,1E-14);

    bool res = g.eval();
    cout << SHOW(res) << endl;
    g.print();


  }
*/

  point3<double> a,b,c;
  a = t.pi[3] - t.pi[0];
  c = t.pi[1] - t.pi[0];
  b = t.pi[2] - t.pi[0];

  crossprod< point3<double> > cross;
  point3<double> N0,N1,N2,N3;
  cross(N0,a-b,c-b);
  cross(N1,a,b);
  cross(N2,a,c);
  cross(N3,b,c);

  point3<double> cc0 = t.circumcenter(0);
  point3<double> cc1 = t.circumcenter(1);
  point3<double> cc2 = t.circumcenter(2);
  point3<double> cc3 = t.circumcenter(3);

  gX.push_back( new gobjglColor3ub(220,20,60) );
  gX.push_back( new gobjMySphereDraw(cc0,gd) );
  gX.push_back( new gobjMySphereDraw(cc1,gd) );
  gX.push_back( new gobjMySphereDraw(cc2,gd) );
  gX.push_back( new gobjMySphereDraw(cc3,gd) );





/*
  // Looking at perpendicular line through midpoints.
  gX.push_back( new gobjglColor3ub(255,0,0) );
  gX.push_back( new gobjglBegin(GL_LINES) );
  gX.push_back( new gobjglVertex3f(m0+N0*-1.0) );
  gX.push_back( new gobjglVertex3f(m0+N0*1.0) );
  gX.push_back( new gobjglVertex3f(m1+N1*-1.0) );
  gX.push_back( new gobjglVertex3f(m1+N1*1.0) );
  gX.push_back( new gobjglVertex3f(m2+N2*-1.0) );
  gX.push_back( new gobjglVertex3f(m2+N2*1.0) );
  gX.push_back( new gobjglVertex3f(m3+N3*-1.0) );
  gX.push_back( new gobjglVertex3f(m3+N3*1.0) );
  gX.push_back( new gobjglEnd() );
*/


/*
  {
    unsigned int n=9;
    double z[(n+1)*n];
    for (unsigned int i=0; i<(n+1)*n; ++i)
      z[i] = 0;

    unsigned int k = 0;

    z[k+0] = a.x;  z[k+1] = c.x;  z[k+2] = -b.x;  z[k+3] = -a.x;  
    z[k+8] = N1.x/N1.dot() - N2.x/N2.dot();
    k += n;

    z[k+0] = a.y;  z[k+1] = c.y;  z[k+2] = -b.y;  z[k+3] = -a.y;  
    z[k+8] = N1.y/N1.dot() - N2.y/N2.dot();
    k += n;

    z[k+0] = a.z;  z[k+1] = c.z;  z[k+2] = -b.z;  z[k+3] = -a.z;  
    z[k+8] = N1.z/N1.dot() - N2.z/N2.dot();
    k += n;


    z[k+4] = a.x - b.x;  z[k+5] = c.x - b.x;  z[k+6] = -b.x;  z[k+7] = c.x;
    z[k+8] = N3.x/N3.dot() - N0.z/N0.dot();
    k += n;

    z[k+4] = a.y - b.y;  z[k+5] = c.y - b.y;  z[k+6] = -b.y;  z[k+7] = c.y;
    z[k+8] = N3.y/N3.dot() - N0.y/N0.dot();
    k += n;

    z[k+4] = a.z - b.z;  z[k+5] = c.z - b.z;  z[k+6] = -b.z;  z[k+7] = c.z;
    z[k+8] = N3.z/N3.dot() - N0.z/N0.dot();
    k += n;


    z[k+0] = a.x;  z[k+1] = c.x;  z[k+6] = -b.x;  z[k+7] = -c.x; 
    z[k+8] = N1.x/N1.dot() - N0.x/N0.dot();
    k += n;

    z[k+0] = a.y;  z[k+1] = c.y;  z[k+6] = -b.y;  z[k+7] = -c.y; 
    z[k+8] = N1.y/N1.dot() - N0.y/N0.dot();
    k += n;

    z[k+0] = a.z;  z[k+1] = c.z;  z[k+6] = -b.z;  z[k+7] = -c.z; 
    z[k+8] = N1.z/N1.dot() - N0.z/N0.dot();

    gausselim<double,10> g(z,1E-14);

    bool res = g.eval();
    cout << SHOW(res) << endl;
    g.print();
  }
*/






}

#endif