//#include <CGAL/Simple_cartesian.h>
//#include <CGAL/Filtered_kernel.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

#include <CGAL/Triangulation_3.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Triangulation_hierarchy_3.h>

#include <CGAL/IO/Geomview_stream.h>
#include <CGAL/IO/Triangulation_geomview_ostream_3.h>

#include <iostream>
#include <fstream>
#include <unistd.h>
#include <vector>

//typedef CGAL::Filtered_kernel<CGAL::Simple_cartesian<double> > K;
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;

typedef K::FT NT;

typedef CGAL::Triangulation_vertex_base_3<K>             Vb;
typedef CGAL::Triangulation_hierarchy_vertex_base_3<Vb>  Vbh;
typedef CGAL::Triangulation_data_structure_3<Vbh>        Tds;
typedef CGAL::Delaunay_triangulation_3<K,Tds>            Dt;
typedef CGAL::Triangulation_hierarchy_3<Dt>              Dh;

typedef Dh::Vertex_iterator Vertex_iterator;
typedef Dh::Vertex_handle   Vertex_handle;
typedef Dh::Cell_handle     Cell_handle;
typedef K::Point_3          Point;

bool is_alpha(const Dh & T,
	      const Dh::Facet & f,
	      NT alpha)
{
  assert( alpha > 0 );
  assert( ! T.is_infinite( f ) );

  Cell_handle c = f.first;
  int ic = f.second;

  Cell_handle n = c->neighbor(ic);

  Point p1, p2, p3;

  if ( ic%2 == 1 ) {
    p1 = c->vertex((ic+1)&3)->point(); // &3 in place of %4
    p2 = c->vertex((ic+2)&3)->point();
    p3 = c->vertex((ic+3)&3)->point();
  }
  else {
    p1 = c->vertex((ic+2)&3)->point();
    p2 = c->vertex((ic+1)&3)->point();
    p3 = c->vertex((ic+3)&3)->point();
  }

  if ( T.is_infinite(n) ) {
    Point Oc = CGAL::circumcenter( c->vertex(0)->point(),
				   c->vertex(1)->point(),
				   c->vertex(2)->point(),
				   c->vertex(3)->point() );
    if ( orientation(p1, p2, p3, Oc) == CGAL::POSITIVE )
      return alpha > CGAL::squared_distance( p1, circumcenter( p1, p2, p3 ) );
    else
      return alpha > CGAL::squared_distance( Oc, p1 );
  }

  if ( T.is_infinite(c) ) {
    Point On = CGAL::circumcenter( n->vertex(0)->point(),
				   n->vertex(1)->point(),
				   n->vertex(2)->point(),
				   n->vertex(3)->point() );
    if ( orientation(p1, p2, p3, On) == CGAL::NEGATIVE )
      return alpha > CGAL::squared_distance( p1, circumcenter( p1, p2, p3 ) );
    else
      return alpha > CGAL::squared_distance( On, p1 );
  }

  // now both c and n are finite

  Point Oc = CGAL::circumcenter( c->vertex(0)->point(),
				 c->vertex(1)->point(),
				 c->vertex(2)->point(),
				 c->vertex(3)->point() );
  Point On = CGAL::circumcenter( n->vertex(0)->point(),
				 n->vertex(1)->point(),
				 n->vertex(2)->point(),
				 n->vertex(3)->point() );

  NT R2c = CGAL::squared_distance( Oc, p1 );
  NT R2n = CGAL::squared_distance( On, p1 );

  if ( orientation( p1,p2,p3,Oc ) != orientation( p1,p2,p3,On ) )
    return alpha < CGAL::max( R2c, R2n ) &&
	   alpha > CGAL::squared_distance( p1, circumcenter( p1,p2,p3 ) );
  else
    return alpha < CGAL::max( R2c, R2n ) &&
           alpha > CGAL::min( R2c, R2n );
}

int main(int argc, char* argv[])
{
  CGAL::Geomview_stream gv(CGAL::Bbox_3(0,0,0, 5, 5, 5));
  gv.set_bg_color(CGAL::Color(0, 200, 200));
  gv.clear();

  Dh T;

  std::cout << "--- Reading file and computing Delaunay... " << std::endl;
  std::ifstream iFile(argv[1], std::ios::in);
  Point p;
  int nbpts = 0;
  while ( iFile >> p ) {
    ++nbpts;
    T.insert( p );
  }

  std::cout << nbpts << " read points" << std::endl;
  std::cout << T.number_of_vertices() << " vertices" << std::endl;

  std::cout << T.number_of_finite_facets()
	    << " facets in Delaunay" << std::endl << std::endl;

  bool again = true;
  double alpha;

  while (again) {
    std::cout << "alpha ? ";
    std::cin >> alpha;
    std::cout << std::endl;

    int nb_facets=0;

    std::cout << "--- computation of the alpha facets and display... " << std::endl;

    std::vector<K::Triangle_3> tris;
    Dh::Finite_facets_iterator tit = T.finite_facets_begin();
    for ( ; tit != T.finite_facets_end() ; ++tit ) {
      if ( is_alpha( T, *tit, alpha ) ) {
	++nb_facets;
	//      gv << T.triangle(*tit);
	tris.push_back(T.triangle(*tit));
      }
    }

    gv.clear();
    gv.draw_triangles(tris.begin(), tris.end());

    std::cout << nb_facets << " alpha facets" << std::endl;

    std::cout << "continue ? yes : 1 / no : 0 ? ";
    std::cin >> again;
    std::cout << std::endl;
  }

  return 0;
}