#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>
int main(int argc, char* argv[]) {
  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 = 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;
}