#include "vtkCGALDelaunay3to2.h"

// VTK related includes
#include "vtkCellArrayIterator.h"
#include "vtkDataSet.h"
#include "vtkDataSetSurfaceFilter.h"
#include "vtkIdList.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkUnstructuredGrid.h"

// CGAL related includes
#include <CGAL/Surface_mesh_default_triangulation_3.h>
#include <CGAL/Surface_mesh.h>
// #include <CGAL/Delaunay_triangulation_3.h>

vtkStandardNewMacro(vtkCGALDelaunay3to2);

using Delaunay     = CGAL::Surface_mesh_default_triangulation_3;
using CGAL_Surface = CGAL::Surface_mesh<Delaunay::Geom_traits::Point_3>;
using Graph_Verts  = boost::graph_traits<CGAL_Surface>::vertex_descriptor;
using Graph_Faces  = boost::graph_traits<CGAL_Surface>::face_descriptor;
using Graph_Coord  = boost::property_map<CGAL_Surface, CGAL::vertex_point_t>::type;

//------------------------------------------------------------------------------
void vtkCGALDelaunay3to2::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);
}

//------------------------------------------------------------------------------
int vtkCGALDelaunay3to2::RequestData(
  vtkInformation*, vtkInformationVector** inputVector, vtkInformationVector* outputVector)
{
  // Get the input and output data objects.
  vtkPolyData* input  = vtkPolyData::GetData(inputVector[0]);
  vtkPolyData* output = vtkPolyData::GetData(outputVector);

  // Create the surface mesh for CGAL
  // --------------------------------

  vtkPoints*    vtkPts     = input->GetPoints();
  vtkDataArray* ptsArr     = vtkPts->GetData();
  const auto    pointRange = vtk::DataArrayTupleRange<3>(ptsArr);

  // CGAL Processing
  // ---------------

  Delaunay delaunay;
  try
  {
    // Mesh points in 3D
    for (auto point : pointRange)
    {
      delaunay.insert(Delaunay::Point(point[0], point[1], point[2]));
    }
  }
  catch (std::exception& e)
  {
    vtkErrorMacro("CGAL Exception: " << e.what());
    return 0;
  }

  // VTK Output
  // ----------

  // Generate 3D mesh:

  vtkNew<vtkPoints> pts;
  const vtkIdType   outNPts = delaunay.number_of_vertices();
  pts->Allocate(outNPts);
  std::map<Delaunay::Point_3, vtkIdType> vmap;
  for (auto vertex : delaunay.all_vertex_handles())
  {
    auto p = vertex->point();
    vtkIdType id =
      pts->InsertNextPoint(CGAL::to_double(p.x()), CGAL::to_double(p.y()), CGAL::to_double(p.z()));
    vmap[p] = id;
  }
  pts->Squeeze();

  // cells
  vtkNew<vtkCellArray> cells;
  cells->AllocateEstimate(delaunay.number_of_facets(), 3);

  for (auto face : delaunay.finite_cell_handles())
  {
    // std::cout << "f: " << face->vertex(0)->point() << std::endl;
    vtkNew<vtkIdList> ids;
    for (int i = 0; i < 4; i++)
    {
      ids->InsertNextId(vmap[face->vertex(i)->point()]);
    }
    cells->InsertNextCell(ids);
  }
  cells->Squeeze();

  vtkNew<vtkUnstructuredGrid> ug;
  ug->SetPoints(pts);
  ug->SetCells(VTK_TETRA, cells);

  // Extract surface in output
  vtkNew<vtkDataSetSurfaceFilter> surface;
  surface->SetInputData(ug);
  surface->Update();

  // output
  output->ShallowCopy(surface->GetOutput());

  return 1;
}
