/**=========================================================================
* @file     main.cpp
*
* @brief    This script reads mesh files and translate them in VTKPolyData.
*           It calculates the main curvature of the mesh, then the shape index in
*           each vertex of the mesh. Generates a Reeb's graph of shape index.
*           Finally, it writes the graph in .vtp format and colors the map with
*           the shape index of mesh.
*
* @authors  Lucie CLERAND, Eve REGA
*           Astrid BEYER
*
* @date     2023
* @version  1.0.1
*
=========================================================================**/

#include <vtkNew.h>
#include <vtkSmartPointer.h>
#include <vtkOBJReader.h>
#include <ttkOFFReader.h>
#include <vtksys/SystemTools.hxx>
#include <vtkPolyData.h>

#include <vtkCurvatures.h>
#include <vtkPointData.h>
#include <vtkCleanPolyData.h>
#include <vtkDoubleArray.h>
#include <vtkIdFilter.h>
#include <ttkFTRGraph.h>
#include <vtkXMLPolyDataWriter.h>
#include <vtkXMLUnstructuredGridWriter.h>

#include <vtkActor.h>
#include <vtkCamera.h>
#include <vtkNamedColors.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>

#include <cassert>
#include <iostream>
#include <cmath>

#include <ttkIcospheresFromPoints.h>
#include <vtkTubeFilter.h>
#include <vtkExtractSurface.h>
#include <ttkGeometrySmoother.h>

#include <vtkDataSetSurfaceFilter.h>
#include <ttkUtils.h>
#include <vtkThresholdPoints.h>
#include <vtkAppendPolyData.h>
#include <vtkSortDataArray.h>

#include "triangulate.h"

int main(int argc, char *argv[])
{
  vtkObject::GlobalWarningDisplayOff();
  if (argc != 2)
  {
    std::cout << "Required arguments: Filename(.obj)" << std::endl;
    return EXIT_FAILURE;
  }

  // Read the input file
  std::string filename = argv[1];
  vtkSmartPointer<vtkPolyData> source;

  std::string extension = vtksys::SystemTools::GetFilenameLastExtension(filename);
  std::transform(extension.begin(), extension.end(), extension.begin(), ::tolower);

  if (extension == ".obj")
  {
    vtkNew<vtkOBJReader> reader;
    reader->SetFileName(filename.c_str());
    reader->Update();
    source = triangulate::buildValidTriangulation(reader->GetOutput());
  }
  else
  {
    std::cout << "Required arguments: Filename(.obj)" << std::endl;
    return EXIT_FAILURE;
  }

  // Compute Max and Min curvature
  vtkSmartPointer<vtkCurvatures> maxCurvaturesFilter = vtkSmartPointer<vtkCurvatures>::New();
  maxCurvaturesFilter->SetInputData(source);
  maxCurvaturesFilter->SetCurvatureTypeToMaximum(); // Set curvature type to maximum (k1)
  maxCurvaturesFilter->Update();
  auto k1 = maxCurvaturesFilter->GetOutput();

  vtkSmartPointer<vtkCurvatures> minCurvaturesFilter = vtkSmartPointer<vtkCurvatures>::New();
  minCurvaturesFilter->SetInputData(source);
  minCurvaturesFilter->SetCurvatureTypeToMinimum(); // Set curvature type to minimum (k2)
  minCurvaturesFilter->Update();
  auto k2 = minCurvaturesFilter->GetOutput();

  // Shape Index computation
  k1->GetPointData()->SetActiveScalars("Maximum_Curvature");
  auto k1Array = k1->GetPointData()->GetAbstractArray("Maximum_Curvature");
  k2->GetPointData()->SetActiveScalars("Minimum_Curvature");
  auto k2Array = k2->GetPointData()->GetAbstractArray("Minimum_Curvature");

  vtkNew<vtkDoubleArray> shapeIndex;
  shapeIndex->SetName("Shape_Index");

  for (vtkIdType i = 0; i < k1->GetNumberOfPoints(); ++i)
  {
    double kmax = k1Array->GetVariantValue(i).ToDouble();
    double kmin = k2Array->GetVariantValue(i).ToDouble();

    if (kmax == kmin)
    {
      shapeIndex->InsertNextTuple1(0.5);
    }
    else
    {
      // double si = (2/M_PI)*atan((kmin+kmax)/(kmin-kmax)); //KOENDERIK ET VAN DOORN [-1,1]
      double si = 0.5 - (1 / M_PI) * atan((kmax + kmin) / (kmin - kmax)); // [0, 1]

      shapeIndex->InsertNextTuple1(si);
    }
  }
  source->GetPointData()->AddArray(shapeIndex);
  source->GetPointData()->SetActiveScalars("Shape_Index");

  /*
    double pourcent = 0.8;

    // Create a vtkThresholdPoints filter for each threshold value
    vtkNew<vtkThresholdPoints> threshold1;
    threshold1->SetInputData(source);
    threshold1->ThresholdBetween(0.0, 0.0 + pourcent*0.25); // Scalar values near 0
    threshold1->Update();

    vtkNew<vtkThresholdPoints> threshold2;
    threshold2->SetInputData(source);
    threshold2->ThresholdBetween(0.5 - pourcent*0.25, 0.5 + pourcent*0.25); // Scalar values near 1
    threshold2->Update();

    vtkNew<vtkThresholdPoints> threshold3;
    threshold3->SetInputData(source);
    threshold3->ThresholdBetween(1.0 - pourcent*0.25, 1.0); // Scalar values near 1
    threshold3->Update();

    // Create a vtkFTRGraph for each thresholded dataset
    vtkNew<ttkFTRGraph> graph1;
    graph1->SetInputData(threshold1->GetOutput());
    graph1->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, "Shape_Index");
    graph1->Update();

    vtkNew<ttkFTRGraph> graph2;
    graph2->SetInputData(threshold2->GetOutput());
    graph2->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, "Shape_Index");
    graph2->Update();

    vtkNew<ttkFTRGraph> graph3;
    graph3->SetInputData(threshold3->GetOutput());
    graph3->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, "Shape_Index");
    graph3->Update();

    // Merge the thresholded datasets together into a single polydata
    vtkNew<vtkAppendPolyData> appendFilter;
    appendFilter->AddInputData(threshold1->GetOutput());
    appendFilter->AddInputData(threshold2->GetOutput());
    appendFilter->AddInputData(threshold3->GetOutput());
    appendFilter->Update();
  */

  // Compute the Reeb Graph of the input dataset
  vtkNew<ttkFTRGraph> reebGraph;
  // reebGraph->SetInputData(appendFilter->GetOutput()); // ttkFTRGraph for the merged polydata
  reebGraph->SetInputData(source);
  reebGraph->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, "Shape_Index");
  reebGraph->Update();

  // Create Icospheres to represent nodes in the Reeb Graph
  vtkNew<ttkIcospheresFromPoints> ttkIcospheresFromPoints{};
  ttkIcospheresFromPoints->SetInputData(0, reebGraph->GetOutput());
  ttkIcospheresFromPoints->SetRadius(0.008);
  ttkIcospheresFromPoints->SetNumberOfSubdivisions(1);
  ttkIcospheresFromPoints->Update();

  // Filters dependencies with Tubes : applies smooth to the geometry of the Reeb Graph
  vtkNew<ttkGeometrySmoother> ttkGeometrySmoother{};
  ttkGeometrySmoother->SetInputConnection(reebGraph->GetOutputPort(1));
  ttkGeometrySmoother->Update();

  vtkNew<vtkGeometryFilter> geometryFilter{};
  geometryFilter->SetInputConnection(ttkGeometrySmoother->GetOutputPort());
  geometryFilter->Update();

  vtkNew<vtkPolyData> polyData{};
  polyData->ShallowCopy(geometryFilter->GetOutput());

  // Create Tubes to represent edges in the Reeb Graph
  vtkNew<vtkTubeFilter> tube{};
  tube->SetInputData(polyData); // smoothed geometry of the Reeb graph
  tube->SetRadius(0.004);
  // tube->SetNumberOfSides(1);
  tube->Update();

  // Save the graph nodes with IcospheresFromPoints
  vtkNew<vtkXMLPolyDataWriter> writerNodes{};
  writerNodes->SetFileName("ReebGraphNodes.vtp");
  writerNodes->SetInputData(ttkIcospheresFromPoints->GetOutput());
  writerNodes->Write();

  // Save the graph edges with Tubes
  vtkNew<vtkXMLPolyDataWriter> writerArcs{};
  writerArcs->SetFileName("ReebGraphArcs.vtp");
  writerArcs->SetInputData(tube->GetOutput());
  writerArcs->Write();

  // Save the graph coloration
  vtkNew<vtkXMLPolyDataWriter> segWriter{};
  segWriter->SetFileName("ShapeIndexMap.vtp");
  segWriter->SetInputConnection(reebGraph->GetOutputPort(2));
  segWriter->Write();

  /*
  // Save the graph nodes (simple graph)
  vtkNew<vtkXMLUnstructuredGridWriter> sWriter{};
  sWriter->SetInputConnection(reebGraph->GetOutputPort(0));
  sWriter->SetFileName("outputNodes.vtp");
  sWriter->Write();

  // Save the graph edges
  vtkNew<vtkXMLUnstructuredGridWriter> sepWriter{};
  sepWriter->SetInputConnection(reebGraph->GetOutputPort(1));
  sepWriter->SetFileName("outputEdges.vtp");
  sepWriter->Write();

  // Save the graph coloration
  vtkNew<vtkXMLPolyDataWriter> segWriter{};
  segWriter->SetInputConnection(reebGraph->GetOutputPort(2));
  segWriter->SetFileName("outputColor.vtp");
  segWriter->Write();*/

  return EXIT_SUCCESS;
}