#define CGAL_MESH_2_OPTIMIZER_VERBOSE
#define CGAL_MESHES_DEBUG_REFINEMENT_POINTS
#define CGAL_MESH_2_DEBUG_BAD_FACES
#define CGAL_MESHES_DEBUG_DOUBLE_MAP
#define CGAL_MESH_2_DEBUG_CONFLICTS_ZONE
#define CGAL_MESH_2_DEBUG_INSERTIONS
#define CGAL_MESH_2_VERBOSE

#include <fstream>

#include "VoronoiGenerator.hpp"

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>

#include <CGAL/Delaunay_mesher_2.h>
#include <CGAL/Delaunay_mesh_face_base_2.h>
#include <CGAL/Delaunay_mesh_vertex_base_2.h>
#include <CGAL/Delaunay_mesh_size_criteria_2.h>

#include <CGAL/Voronoi_diagram_2.h>
#include <CGAL/Delaunay_triangulation_adaptation_traits_2.h>
#include <CGAL/Delaunay_triangulation_adaptation_policies_2.h>

#include <CGAL/Boolean_set_operations_2.h>

#include <CGAL/Cartesian_converter.h>

using namespace dem;

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Triangulation_vertex_base_2<K> Vb;
typedef CGAL::Delaunay_mesh_face_base_2<K> Fb;
typedef CGAL::Triangulation_data_structure_2<Vb, Fb> Tds;
typedef CGAL::Constrained_Delaunay_triangulation_2<K, Tds> CDT;
typedef CGAL::Delaunay_mesh_size_criteria_2<CDT> Criteria;
typedef CGAL::Delaunay_mesher_2<CDT, Criteria> Mesher;

typedef CGAL::Delaunay_triangulation_adaptation_traits_2<CDT>                 AT;
typedef CGAL::Delaunay_triangulation_caching_degeneracy_removal_policy_2<CDT> AP;
typedef CGAL::Voronoi_diagram_2<CDT,AT,AP>                                    VD;

typedef CGAL::Exact_predicates_exact_constructions_kernel exactK;
typedef K::Point_2 Point;
typedef CGAL::Polygon_2<exactK> Polygon;
typedef CGAL::Vector_2<K> Vector_2;
typedef CGAL::Polygon_with_holes_2<exactK> Polygon_with_holes_2;
typedef CGAL::Cartesian_converter<exactK,K> EK_to_IK;


VoronoiGenerator::VoronoiGenerator(std::vector<Point> const& boundaryPtsIn)
  : boundaryPts(boundaryPtsIn)
{

  for(auto& pt : boundaryPts)
    bdryPoly.push_back( exactK::Point_2(pt.x(),pt.y()));

}

VoronoiGenerator::VoronoiGenerator(Polygon const& polygonIn)
  : bdryPoly(polygonIn)
{

  EK_to_IK to_inexact;

  for (auto vi=bdryPoly.vertices_begin(); vi!=bdryPoly.vertices_end(); ++vi)
    boundaryPts.push_back(to_inexact(*vi));

}


std::vector<Polygon> VoronoiGenerator::Generate()
{

  // The constrained Delaunay triangulation (dual of the Voronoi)
  CDT cdt;

  cdt.insert_constraint(boundaryPts.begin(), boundaryPts.end(), true);

  Mesher mesher(cdt);
  mesher.init(false);
  mesher.set_criteria(Criteria(0.25, 1.0));
  std::cout <<"\nNumber of faces: " << cdt.number_of_faces() << "\n" << std::endl;

  std::cout << "\nFirst refinement step:" << std::endl;
  mesher.try_one_step_refine_mesh();
  std::cout <<"\nNumber of faces: " << cdt.number_of_faces() << "\n" << std::endl;

  std::cout << "\nSecond refinement step:" << std::endl;
  std::cout << "Is refinement done: " << mesher.is_refinement_done() << std::endl;
  mesher.try_one_step_refine_mesh();
  std::cout <<"\nNumber of faces: " << cdt.number_of_faces() << "\n" << std::endl;
  /*
  std::cout << "\nThird refinement step:" << std::endl;
  std::cout << "Is refinement done: " << mesher.is_refinement_done() << std::endl;
  mesher.try_one_step_refine_mesh();
  std::cout <<"\nNumber of faces: " << cdt.number_of_faces() << "\n" << std::endl;

  mesher.refine_mesh();
  */
//  CGAL::refine_Delaunay_mesh_2(cdt, Criteria(0.125, 0.05));


  std::vector<Polygon> allCells;
  std::vector<exactK::Point_2> cellVerts;


  for (auto face=cdt.finite_faces_begin(); face!=cdt.finite_faces_end(); ++face)
  {

    std::cout << "Area: " << cdt.triangle(face).area() << std::endl;

    cellVerts.clear();

    cellVerts.push_back(exactK::Point_2(cdt.triangle(face)[0].x(),
					cdt.triangle(face)[0].y()));
    cellVerts.push_back(exactK::Point_2(cdt.triangle(face)[1].x(),
					cdt.triangle(face)[1].y()));
    cellVerts.push_back(exactK::Point_2(cdt.triangle(face)[2].x(),
					cdt.triangle(face)[2].y()));

    Polygon voroCell(cellVerts.begin(), cellVerts.end());
    allCells.push_back(voroCell);

  }


  /*
  // Construct the Voronoi diagram from the Delaunay triangulation
  VD voro(cdt);

  int faceNum = 0;
  std::vector<Polygon> allCells;

  for(auto it = voro.bounded_faces_begin(); it != voro.bounded_faces_end(); ++it, ++faceNum)
  {

    auto edgeIterStart = it->ccb();
    auto edgeIter = edgeIterStart;

    //Polygon_2 voroCell;
    std::vector<exactK::Point_2> cellVerts;

    do {

      Point newPt;

      if(edgeIter->is_segment()){

	std::cout << "a" << std::endl;
	//	std::cout << "Is halfedge:" << it->is_halfedge_on_ccb(*e
	std::cout << "Valid: " << edgeIter->is_unbounded() << std::endl;
	std::cout << "Target: " << edgeIter->target()->point() << std::endl;
	std::cout << "Source: " << edgeIter->has_source() << std::endl;;
	newPt = edgeIter->source()->point();
	cellVerts.push_back(exactK::Point_2(newPt.x(), newPt.y()));

	newPt = edgeIter->target()->point();
	cellVerts.push_back(exactK::Point_2(newPt.x(), newPt.y()));

      }else if (edgeIter->is_ray()){

	std::cout << "b" << std::endl;

	Point p1 = edgeIter->face()->dual()->point();
	Point p2 = edgeIter->twin()->face()->dual()->point();
	Point halfPt = p1 + 0.5*(p2-p1);

	if(edgeIter->has_target()){

	  std::cout << "c" << std::endl;
	  Point tgtPt = edgeIter->target()->point();
	  newPt = tgtPt + bigNum*(halfPt - tgtPt);

	  cellVerts.push_back(exactK::Point_2(newPt.x(), newPt.y()));
	  cellVerts.push_back(exactK::Point_2(tgtPt.x(), tgtPt.y()));

	}else if(edgeIter->has_source()){

	  std::cout << "d" << std::endl;
	  Point srcPt = edgeIter->source()->point();
	  newPt = srcPt + bigNum*(halfPt - srcPt);

	  cellVerts.push_back(exactK::Point_2(srcPt.x(), srcPt.y()));
	  cellVerts.push_back(exactK::Point_2(newPt.x(), newPt.y()));

	}
      }

    } while ( ++edgeIter != edgeIterStart );


    cellVerts.erase( std::unique(cellVerts.begin(), cellVerts.end()), cellVerts.end());

    if(*(cellVerts.end()-1) == cellVerts.at(0))
      cellVerts.pop_back();

    // Copy the vertices into a polygon
    Polygon voroCell(cellVerts.begin(), cellVerts.end());

    // Get the intersection of the
    std::vector<Polygon_with_holes_2> intersection;
    CGAL::intersection (voroCell, bdryPoly, std::back_inserter(intersection));

    assert(intersection.size()==1);
    allCells.push_back(intersection.at(0).outer_boundary());

    allCells.push_back(voroCell);


  } // end of loop over Voronoi cells
*/
  return allCells;

}


void dem::VoronoiToFile(std::vector<Polygon>& polys, std::string filename)
{

  std::ofstream polysOut(filename);

  if (polysOut.is_open())
  {

    for (auto p : polys)
    {

      polysOut << p.size() << std::endl;

      for (auto it=p.vertices_begin(); it!=p.vertices_end(); ++it)
	polysOut << it->x() << "  " << it->y() << std::endl;

    }

  }

}