/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkLSCMFilter.h

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/
/**
 * @class   vtkLSCMFilter
 * @brief   Compute least-squares comformal maps using chart ids on triangles.
 *
 * This filter takes input as a triangular surface mesh with chart assignment
 * for each triangle and outputs: (1) a vtkPartitionedDataSet of vtkPolyData
 * objects holding triangles of each chart with u-v coordinates and input
 * global id assigned to each point, and (2) a vtkPartitionedDataSet of
 * vtkPolyData objects holding polyline cells (VTK_POLY_LINE) corresponding to
 * seams of the chart and boundaries of the polydata. Each polydata have two
 * vtkIdTypeArray as cell data that mark each polyline with partition ID and
 * cell ID of its partner on another output polydata. Polylines representing
 * input boundaries have no partner and will have -1 for both cell data.
 *
 * The computation of LSCM is a least square problem with linear constraints
 * imposed by two anchors (points that are farthest in the input polydata
 * in the same chart). The problem is formulated into KKT equations and is
 * solved using QR decomposition with column pivoting. Within each chart, two
 * points that are farthest from each other in the input polydata are used as
 * the anchor points with anchorPoint0 assigned to origin and anchorPoint1
 * assigned to (0, originalDistance)
 *
 * The current computation of LSCM is implemented in sequential loops with
 * dense matrix operations using Eigen library.
 */

#ifndef vtk_aeva_ext_vtkLSCMFilter_h
#define vtk_aeva_ext_vtkLSCMFilter_h

#include "vtk/aeva/ext/AEVAExtModule.h"
#include "vtkPartitionedDataSet.h"
#include "vtkPolyDataAlgorithm.h"
#include "vtkeigen/eigen/Dense"

using namespace std;
using namespace Eigen;

class AEVAEXT_EXPORT vtkLSCMFilter : public vtkPolyDataAlgorithm
{
public:
  vtkTypeMacro(vtkLSCMFilter, vtkPolyDataAlgorithm);
  void PrintSelf(ostream& os, vtkIndent indent) override;

  vtkLSCMFilter(const vtkLSCMFilter&) = delete;
  vtkLSCMFilter& operator=(const vtkLSCMFilter&) = delete;

  /**
   * Construct object with
   */
  static vtkLSCMFilter* New();

  /// Used to index input-data arrays
  enum InputArrays
  {
    CHART_ID = 0 //!< The index of the input array holding chart IDs (must be cell-data).
  };

  /// Used to index output datasets
  enum OutputPorts
  {
    ATLAS = 0,    //!< The output texture atlas.
    BOUNDARY = 1, //!< Shared output chart boundaries.
    SEGMENTS = 2  //!< Original individual segments based on chart ids.
  };

  /// All outputs are partitioned data; provide methods that return this type:
  vtkPartitionedDataSet* GetOutput(int port = 0);
  vtkPartitionedDataSet* GetAtlasOutput() { return this->GetOutput(ATLAS); }
  vtkPartitionedDataSet* GetBoundaryOutput() { return this->GetOutput(BOUNDARY); }
  vtkPartitionedDataSet* GetSegmentsOutput() { return this->GetOutput(SEGMENTS); }

private:
  //@{
  /**
   * Compute cotangent weighted discrete Laplacian-Beltrami operator
   */
  static void ComputeCotMatrix(const vtkSmartPointer<vtkPolyData>& inputPD, MatrixXd& lapMat);
  //@}

  //@{
  /**
   * Compute Area matrix
   */
  static void ComputeAreaMatrix(const vtkSmartPointer<vtkPolyData>& inputPD, MatrixXd& AMat2);
  //@}

  /**
   * Compute LSCM solution
   */
  static void ComputeLSCM(VectorXd& lVec2,
    MatrixXd& lapMat2,
    MatrixXd& aMat2,
    const vtkSmartPointer<vtkIdList>& anchorIds,
    VectorXd& anchorDisps);
  //@}

  /*
   * Functor for computing local LSCM coordinates in parallel
   */
  class ComputeLocalCoordinates
  {
  public:
    explicit ComputeLocalCoordinates(vtkSmartPointer<vtkPartitionedDataSet> outputPDS0,
      vtkSmartPointer<vtkPartitionedDataSet> outputPDS1,
      vector<double> chartAreas)
      : m_outputPDS0(std::move(outputPDS0))
      , m_outputPDS1(std::move(outputPDS1))
      , m_chartAreas(std::move(chartAreas))
    {
    }
    vtkSmartPointer<vtkPartitionedDataSet> m_outputPDS0;
    vtkSmartPointer<vtkPartitionedDataSet> m_outputPDS1;
    vector<double> m_chartAreas;
    void operator()(vtkIdType begin, vtkIdType end);
  };

protected:
  vtkLSCMFilter();
  ~vtkLSCMFilter() override = default;

  int RequestData(vtkInformation* request,
    vtkInformationVector** inputVector,
    vtkInformationVector* outputVector) override;
  int FillOutputPortInformation(int port, vtkInformation* info) override;
};
#endif
