#include <iostream>
using namespace std;

#include <delaunay3D.h>
#include <halfspaceD3.h>
#include <tetrahedron.h>


delaunay3D::delaunay3D
( 
  vector<point3<double> > const & vi_ 
)
  : vi(vi_)
{
  n = vi.size();
}

void delaunay3D::eval()
{
  tet.clear();

  if (n>200)
  {
    cout << "error:  n is too large." << endl;
    cout << "  this is an O(n^5) algorithm." << endl;
    return;
  }

  uint a,i,j,k;

  for (i=0; i<n-3; ++i)
  {
    for (j=i+1; j<n-2; ++j)
    {
      for (k=j+1; k<n-1; ++k)
      {
        for (a=k+1; a<n; ++a)
        {
          evaluateTetrahedron(i,j,k,a);
        }
      }
    }
  }
}

bool const delaunay3D::evaluateTetrahedron
(
  uintc i,
  uintc j,
  uintc k,
  uintc a
)
{
  point3<double> const & P0(vi[i]);
  point3<double> const & P1(vi[j]);
  point3<double> const & P2(vi[k]);
  point3<double> const & P3(vi[a]);

  //tetrahedron< point3<double>,double> const t(P0,P1,P2,P3);
  tetrahedron< point3<double>,double> t;
  t.constructUnordered(P0,P1,P2,P3);
  
  point3<double> c0;
  t.circumcenter(c0);
  point3<double> const c(c0);

  // Avoid the square root function by comparing 
  //   the distance squared.
  double const dist = (P0-c).dot();

  for (uint w=0; w<n; ++w)
  {
    if (w==i)
      continue;
    if (w==j)
      continue;
    if (w==k)
      continue;
    if (w==a)
      continue;

    if ( ((vi[w]-c).dot()) < dist )
      return false;
  }

  halfspaceD3< point3<double>,double> h(P0,P1,P2);

  //if (! h.eval(P3))
  if (! h.isInside(P3))
    tet.push_back( point4<uint>(i,j,k,a) );
  else
    tet.push_back( point4<uint>(k,j,i,a) );

  return true;
}