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

  Program:   Visualization Toolkit
  Module:    vtkReebGraph.cxx

  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.

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

/*----------------------------------------------------------------------------
 Copyright (c) Sandia Corporation
 See Copyright.txt or http://www.paraview.org/HTML/Copyright.html for details.
----------------------------------------------------------------------------*/

#include "vtkReebGraph.h"

#include "vtkCell.h"
#include "vtkDataArray.h"
#include "vtkEdgeListIterator.h"
#include "vtkIdTypeArray.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkReebGraphSimplificationMetric.h"
#include "vtkUnstructuredGrid.h"
#include "vtkVariantArray.h"

#include <algorithm>
#include <map>
#include <queue>
#include <vector>

//------------------------------------------------------------------------------
// Contain all of the internal data structures, and macros, in the
// implementation.
namespace
{
//------------------------------------------------------------------------------
inline bool vtkReebGraphVertexSoS(
  const std::pair<int, double>& v0, const std::pair<int, double>& v1)
{
  return ((v0.second < v1.second) || ((v0.second == v1.second) && (v0.first < v1.first)));
}
}

// INTERNAL MACROS ---------------------------------------------------------
#define vtkReebGraphSwapVars(type, var1, var2)                                                     \
  {                                                                                                \
    type tmp;                                                                                      \
    tmp = (var1);                                                                                  \
    (var1) = (var2);                                                                               \
    (var2) = tmp;                                                                                  \
  }

#define vtkReebGraphInitialStreamSize 1000

#define vtkReebGraphIsSmaller(myReebGraph, nodeId0, nodeId1, node0, node1)                         \
  ((node0->Value < node1->Value) || (node0->Value == node1->Value && (nodeId0) < (nodeId1)))

#define vtkReebGraphGetArcPersistence(rg, a)                                                       \
  ((this->GetNode(a->NodeId1)->Value - this->GetNode(a->NodeId0)->Value) /                         \
    (this->MaximumScalarValue - this->MinimumScalarValue))

#define vtkReebGraphIsHigherThan(rg, N0, N1, n0, n1)                                               \
  ((n0->Value > n1->Value) || (n0->Value == n1->Value && n0->VertexId > n1->VertexId))

// Note: usually this macro is called after the node has been finalized,
// otherwise the behaviour is undefined.
#define vtkReebGraphIsRegular(rg, n)                                                               \
  ((!(n)->IsCritical) &&                                                                           \
    ((n)->ArcDownId && !this->GetArc((n)->ArcDownId)->ArcDwId1 && (n)->ArcUpId &&                  \
      !this->GetArc((n)->ArcUpId)->ArcDwId0))

#define vtkReebGraphAddUpArc(rg, N, A)                                                             \
  {                                                                                                \
    vtkReebNode* n = this->GetNode(N);                                                             \
    vtkReebArc* a = this->GetArc(A);                                                               \
    a->ArcUpId0 = 0;                                                                               \
    a->ArcDwId0 = n->ArcUpId;                                                                      \
    if (n->ArcUpId)                                                                                \
      this->GetArc(n->ArcUpId)->ArcUpId0 = (A);                                                    \
    n->ArcUpId = (A);                                                                              \
  }

#define vtkReebGraphAddDownArc(rg, N, A)                                                           \
  {                                                                                                \
    vtkReebNode* n = this->GetNode(N);                                                             \
    vtkReebArc* a = this->GetArc(A);                                                               \
    a->ArcUpId1 = 0;                                                                               \
    a->ArcDwId1 = n->ArcDownId;                                                                    \
    if (n->ArcDownId)                                                                              \
      this->GetArc(n->ArcDownId)->ArcUpId1 = (A);                                                  \
    n->ArcDownId = (A);                                                                            \
  }

#define vtkReebGraphRemoveUpArc(rg, N, A)                                                          \
  {                                                                                                \
    vtkReebNode* n = this->GetNode(N);                                                             \
    vtkReebArc* a = this->GetArc(A);                                                               \
    if (a->ArcUpId0)                                                                               \
      this->GetArc(a->ArcUpId0)->ArcDwId0 = a->ArcDwId0;                                           \
    else                                                                                           \
      n->ArcUpId = a->ArcDwId0;                                                                    \
    if (a->ArcDwId0)                                                                               \
      this->GetArc(a->ArcDwId0)->ArcUpId0 = a->ArcUpId0;                                           \
  }

#define vtkReebGraphRemoveDownArc(rg, N, A)                                                        \
  {                                                                                                \
    vtkReebNode* n = this->GetNode(N);                                                             \
    vtkReebArc* a = this->GetArc(A);                                                               \
    if (a->ArcUpId1)                                                                               \
      this->GetArc(a->ArcUpId1)->ArcDwId1 = a->ArcDwId1;                                           \
    else                                                                                           \
      n->ArcDownId = a->ArcDwId1;                                                                  \
    if (a->ArcDwId1)                                                                               \
      this->GetArc(a->ArcDwId1)->ArcUpId1 = a->ArcUpId1;                                           \
  }

#ifndef vtkReebGraphMax
#define vtkReebGraphMax(a, b) (((a) >= (b)) ? (a) : (b))
#endif

#define vtkReebGraphStackPush(N)                                                                   \
  {                                                                                                \
    if (nstack == mstack)                                                                          \
    {                                                                                              \
      mstack = vtkReebGraphMax(128, mstack * 2);                                                   \
      int* oldstack = stack;                                                                       \
      stack = (int*)realloc(stack, sizeof(int) * mstack);                                          \
      if (!stack)                                                                                  \
      {                                                                                            \
        free(oldstack);                                                                            \
        assert(0 && "Ran out of memory");                                                          \
      }                                                                                            \
    }                                                                                              \
    stack[nstack++] = (N);                                                                         \
  }

#define vtkReebGraphStackSize() (nstack)

#define vtkReebGraphStackTop() (stack[nstack - 1])

#define vtkReebGraphStackPop() (--nstack)

//------------------------------------------------------------------------------
// PIMPLed classes...
class vtkReebGraph::Implementation
{
public:
  Implementation()
  {
    this->historyOn = false;

    this->MainNodeTable.Buffer = (vtkReebNode*)malloc(sizeof(vtkReebNode) * 2);
    this->MainArcTable.Buffer = (vtkReebArc*)malloc(sizeof(vtkReebArc) * 2);
    this->MainLabelTable.Buffer = (vtkReebLabel*)malloc(sizeof(vtkReebLabel) * 2);

    this->MainNodeTable.Size = 2;
    this->MainNodeTable.Number = 1; // the item "0" is blocked
    this->MainArcTable.Size = 2;
    this->MainArcTable.Number = 1;
    this->MainLabelTable.Size = 2;
    this->MainLabelTable.Number = 1;

    this->MainNodeTable.FreeZone = 1;
    // Clear node
    this->GetNode(1)->ArcUpId = ((int)-2);
    // Initialize DownArc
    this->GetNode(1)->ArcDownId = 0;
    this->MainArcTable.FreeZone = 1;
    // Clear Arc label 1
    this->GetArc(1)->LabelId1 = ((int)-2);
    // Initialize Arc label 0
    this->GetArc(1)->LabelId0 = 0;
    this->MainLabelTable.FreeZone = 1;
    // Clear label
    this->GetLabel(1)->HNext = ((int)-2);
    // Initialize Arc id
    this->GetLabel(1)->ArcId = 0;

    this->MinimumScalarValue = 0;
    this->MaximumScalarValue = 0;

    this->ArcNumber = 0;
    this->NodeNumber = 0;
    this->LoopNumber = 0;
    this->RemovedLoopNumber = 0;
    this->ArcLoopTable = nullptr;

    this->currentNodeId = 0;
    this->currentArcId = 0;

    // streaming support
    this->VertexMapSize = 0;
    this->VertexMapAllocatedSize = 0;
    this->TriangleVertexMapSize = 0;
    this->TriangleVertexMapAllocatedSize = 0;
  }
  ~Implementation()
  {
    free(this->MainNodeTable.Buffer);
    this->MainNodeTable.Buffer = nullptr;

    free(this->MainArcTable.Buffer);
    this->MainArcTable.Buffer = nullptr;

    free(this->MainLabelTable.Buffer);
    this->MainLabelTable.Buffer = nullptr;

    this->MainNodeTable.Size = this->MainNodeTable.Number = 0;
    this->MainArcTable.Size = this->MainArcTable.Number = 0;
    this->MainLabelTable.Size = this->MainLabelTable.Number = 0;

    this->MainNodeTable.FreeZone = 0;
    this->MainArcTable.FreeZone = 0;
    this->MainLabelTable.FreeZone = 0;

    if (this->ArcLoopTable)
      free(this->ArcLoopTable);

    if (this->VertexMapAllocatedSize)
      free(this->VertexMap);

    if (this->TriangleVertexMapAllocatedSize)
      free(this->TriangleVertexMap);
  }

  typedef unsigned long long vtkReebLabelTag;

  struct vtkReebCancellation_t
  {
    std::vector<std::pair<int, int>> removedArcs;
    std::vector<std::pair<int, int>> insertedArcs;
  };
  using vtkReebCancellation = struct vtkReebCancellation_t;

  // Node structure
  struct vtkReebNode_t
  {
    vtkIdType VertexId;
    double Value;
    vtkIdType ArcDownId;
    vtkIdType ArcUpId;
    bool IsFinalized;
    bool IsCritical;
  };
  using vtkReebNode = struct vtkReebNode_t;

  // Arc structure
  struct vtkReebArc_t
  {
    vtkIdType NodeId0, ArcUpId0, ArcDwId0;
    vtkIdType NodeId1, ArcUpId1, ArcDwId1;
    vtkIdType LabelId0, LabelId1;
  };
  using vtkReebArc = struct vtkReebArc_t;

  // Label structure
  struct vtkReebLabel_t
  {
    vtkIdType ArcId;
    vtkIdType HPrev, HNext; // "horizontal" (for a single arc)
    vtkReebLabelTag label;
    vtkIdType VPrev, VNext; // "vertical" (for a sequence of arcs)
  };
  using vtkReebLabel = struct vtkReebLabel_t;

  struct vtkReebPath
  {
    double SimplificationValue;
    int ArcNumber;
    vtkIdType* ArcTable;
    int NodeNumber;
    vtkIdType* NodeTable;

    inline bool operator<(struct vtkReebPath const& E) const
    {
      return !((SimplificationValue < E.SimplificationValue) ||
        (SimplificationValue == E.SimplificationValue && ArcNumber < E.ArcNumber) ||
        (SimplificationValue == E.SimplificationValue && ArcNumber == E.ArcNumber &&
          NodeTable[NodeNumber - 1] < E.NodeTable[E.NodeNumber - 1]));
      /*      return !((
              (MaximumScalarValue - MinimumScalarValue)
                < (E.MaximumScalarValue - E.MinimumScalarValue)) ||
                 ((MaximumScalarValue - MinimumScalarValue)
                   == (E.MaximumScalarValue-E.MinimumScalarValue)
                     && ArcNumber < E.ArcNumber) ||
                 ((MaximumScalarValue - MinimumScalarValue)
                   == (E.MaximumScalarValue - E.MinimumScalarValue)
                     && ArcNumber == E.ArcNumber
                       && NodeTable[NodeNumber - 1]<E.NodeTable[E.NodeNumber - 1])
               );*/
    }
  };

  struct
  {
    int Size, Number, FreeZone;
    vtkReebArc* Buffer;
  } MainArcTable;

  struct
  {
    int Size, Number, FreeZone;
    vtkReebNode* Buffer;
  } MainNodeTable;

  struct
  {
    int Size, Number, FreeZone;
    vtkReebLabel* Buffer;
  } MainLabelTable;

  vtkReebPath FindPath(
    vtkIdType arcId, double simplificationThreshold, vtkReebGraphSimplificationMetric* metric);

  // INTERNAL METHODS --------------------------------------------------------

  vtkIdType AddArc(vtkIdType nodeId0, vtkIdType nodeId1);

  // Description:
  // Get the node specified from the graph.
  vtkReebNode* GetNode(vtkIdType nodeId);

  // Description:
  // Get the arc specified from the graph.
  vtkReebArc* GetArc(vtkIdType arcId);

  // Description:
  // Get the Label specified from the graph.
  vtkReebLabel* GetLabel(vtkIdType labelId);

  // Description:
  // Collapse (consolidate) two nodes in the graph.
  void CollapseVertex(vtkIdType N, vtkReebNode* n);

  // Description:
  // Triggers customized code for simplification metric evaluation.
  double ComputeCustomMetric(vtkReebGraphSimplificationMetric* simplificationMetric, vtkReebArc* a);

  // Description:
  //  Add a monotonic path between nodes.
  vtkIdType AddPath(int nodeNumber, vtkIdType* nodeOffset, vtkReebLabelTag label);

  // Description:
  //   Add a vertex from the mesh to the Reeb graph.
  vtkIdType AddMeshVertex(vtkIdType vertexId, double scalar);

  // Description:
  //   Add a triangle from the mesh to the Reeb grpah.
  int AddMeshTriangle(
    vtkIdType vertex0Id, double f0, vtkIdType vertex1Id, double f1, vtkIdType vertex2Id, double f2);

  // Description:
  //   Add a tetrahedron from the mesh to the Reeb grpah.
  int AddMeshTetrahedron(vtkIdType vertex0Id, double f0, vtkIdType vertex1Id, double f1,
    vtkIdType vertex2Id, double f2, vtkIdType vertex3Id, double f3);

  // Description:
  // "Zip" the corresponding paths when the interior of a simplex is added to
  // the Reeb graph.
  void Collapse(vtkIdType startingNode, vtkIdType endingNode, vtkReebLabelTag startingLabel,
    vtkReebLabelTag endingLabel);

  // Description:
  // Finalize a vertex.
  void EndVertex(vtkIdType N);

  // Description:
  // Remove an arc during filtering by persistence.
  void FastArcSimplify(vtkIdType arcId, int arcNumber, vtkIdType* arcTable);

  // Description:
  // Remove arcs below the provided persistence.
  int SimplifyBranches(
    double simplificationThreshold, vtkReebGraphSimplificationMetric* simplificationMetric);

  // Description:
  // Remove the loops below the provided persistence.
  int SimplifyLoops(
    double simplificationThreshold, vtkReebGraphSimplificationMetric* simplificationMetric);

  // Description:
  // Update the vtkMutableDirectedGraph internal structure after filtering, with
  // deg-2 nodes maintaining.
  int CommitSimplification();

  // Description:
  // Retrieve downwards labels.
  vtkIdType FindDwLabel(vtkIdType nodeId, vtkReebLabelTag label);

  // Description
  // Find greater arc (persistence-based simplification).
  vtkIdType FindGreater(vtkIdType nodeId, vtkIdType startingNodeId, vtkReebLabelTag label);

  // Description:
  // Find corresponding joining saddle node (persistence-based simplification).
  vtkIdType FindJoinNode(vtkIdType arcId, vtkReebLabelTag label, bool onePathOnly = false);

  // Description:
  // Find smaller arc (persistence-based simplification).
  vtkIdType FindLess(vtkIdType nodeId, vtkIdType startingNodeId, vtkReebLabelTag label);

  // Description:
  // Compute the loops in the Reeb graph.
  void FindLoops();

  // Description:
  // Find corresponding splitting saddle node (persistence-based
  // simplification).
  vtkIdType FindSplitNode(vtkIdType arcId, vtkReebLabelTag label, bool onePathOnly = false);

  // Description:
  // Retrieve upwards labels.
  vtkIdType FindUpLabel(vtkIdType nodeId, vtkReebLabelTag label);

  // Description:
  // Flush labels.
  void FlushLabels();

  // Description:
  // Resize the arc table.
  void ResizeMainArcTable(int newSize);

  // Description:
  // Resize the label table.
  void ResizeMainLabelTable(int newSize);

  // Description:
  // Resize the node table.
  void ResizeMainNodeTable(int newSize);

  // Description:
  // Set a label.
  void SetLabel(vtkIdType A, vtkReebLabelTag Label);

  // Description:
  // Simplify labels.
  void SimplifyLabels(
    vtkIdType nodeId, vtkReebLabelTag onlyLabel = 0, bool goDown = true, bool goUp = true);

  // ACCESSORS

  // Description:
  // Returns the Id of the lower node of the arc specified by 'arcId'.
  vtkIdType GetArcDownNodeId(vtkIdType arcId);

  // Description:
  // Return the Id of the upper node of the arc specified by 'arcId'.
  vtkIdType GetArcUpNodeId(vtkIdType arcId);

  // Description:
  // Iterates forwards through the arcs of the Reeb graph.
  //
  // The first time this method is called, the first arc's Id will be returned.
  // When the last arc is reached, this method will keep on returning its Id at
  // each call. See 'GetPreviousArcId' to go back in the list.
  vtkIdType GetNextArcId();

  // Description:
  // Iterates forwards through the nodes of the Reeb graph.
  //
  // The first time this method is called, the first node's Id will be returned.
  // When the last node is reached, this method will keep on returning its Id at
  // each call. See 'GetPreviousNodeId' to go back in the list.
  vtkIdType GetNextNodeId();

  // Description:
  // Copy into 'arcIdList' the list of the down arcs' Ids, given a node
  // specified by 'nodeId'.
  void GetNodeDownArcIds(vtkIdType nodeId, vtkIdList* arcIdList);

  // Description:
  // Returns the scalar field value of the node specified by 'nodeId'.
  double GetNodeScalarValue(vtkIdType nodeId);

  // Description:
  // Copy into 'arcIdList' the list of the up arcs' Ids, given a node specified
  // by 'nodeId'.
  void GetNodeUpArcIds(vtkIdType nodeId, vtkIdList* arcIdList);

  // Description:
  // Returns the corresponding vertex Id (in the simplicial mesh, vtkPolyData),
  // given a node specified by 'nodeId'.
  vtkIdType GetNodeVertexId(vtkIdType nodeId);

  // Description:
  // Returns the number of arcs in the Reeb graph.
  int GetNumberOfArcs();

  // Description:
  // Returns the number of connected components of the Reeb graph.
  int GetNumberOfConnectedComponents();

  // Description:
  // Returns the number of nodes in the Reeb graph.
  int GetNumberOfNodes();

  // Description:
  // Returns the number of loops (cycles) in the Reeb graph.
  //
  // Notice that for closed PL 2-manifolds, this number equals the genus of the
  // manifold.
  //
  // Reference:
  // "Loops in Reeb graphs of 2-manifolds",
  // K. Cole-McLaughlin, H. Edelsbrunner, J. Harer, V. Natarajan, and V.
  // Pascucci,
  // ACM Symposium on Computational Geometry, pp. 344-350, 2003.
  int GetNumberOfLoops();

  // Description:
  // Iterates backwards through the arcs of the Reeb graph.
  //
  // When the first arc is reached, this method will keep on returning its Id at
  // each call. See 'GetNextArcId' to go forwards in the list.
  vtkIdType GetPreviousArcId();

  // Description:
  // Iterates backwards through the nodes of the Reeb graph.
  //
  // When the first node is reached, this method will keep on returning its Id
  // at each call. See 'GetNextNodeId' to go forwards in the list.
  vtkIdType GetPreviousNodeId();

  // Description:
  // Implementations of the stream classes to operate on the streams...
  int StreamTetrahedron(vtkIdType vertex0Id, double scalar0, vtkIdType vertex1Id, double scalar1,
    vtkIdType vertex2Id, double scalar2, vtkIdType vertex3Id, double scalar3);
  int StreamTriangle(vtkIdType vertex0Id, double scalar0, vtkIdType vertex1Id, double scalar1,
    vtkIdType vertex2Id, double scalar2);

  void DeepCopy(Implementation* src);

  // Description:
  // Data storage.
  std::map<int, int> VertexStream;
  std::vector<vtkReebCancellation> cancellationHistory;
  // Streaming support
  int VertexMapSize, VertexMapAllocatedSize, TriangleVertexMapSize, TriangleVertexMapAllocatedSize;
  bool historyOn;

  vtkIdType* VertexMap;
  int* TriangleVertexMap;

  double MinimumScalarValue, MaximumScalarValue;

  // Arcs and nodes
  int ArcNumber, NodeNumber;

  // Loops
  int LoopNumber, RemovedLoopNumber;
  vtkIdType* ArcLoopTable;

  // CC
  int ConnectedComponentNumber;

  std::map<int, double> ScalarField;

  vtkIdType currentNodeId, currentArcId;

  vtkDataSet* inputMesh;
  vtkDataArray* inputScalarField;

  vtkReebGraph* Parent;
};

//------------------------------------------------------------------------------
vtkReebGraph::Implementation::vtkReebNode* vtkReebGraph::Implementation::GetNode(vtkIdType nodeId)
{
  return (this->MainNodeTable.Buffer + nodeId);
}

//------------------------------------------------------------------------------
vtkReebGraph::Implementation::vtkReebArc* vtkReebGraph::Implementation::GetArc(vtkIdType arcId)
{
  return (this->MainArcTable.Buffer + arcId);
}

//------------------------------------------------------------------------------
vtkReebGraph::Implementation::vtkReebLabel* vtkReebGraph::Implementation::GetLabel(
  vtkIdType labelId)
{
  return (this->MainLabelTable.Buffer + labelId);
}

//------------------------------------------------------------------------------
void vtkReebGraph::Implementation::CollapseVertex(vtkIdType N, vtkReebNode* n)
{
  int Lb, Lnext, La;
  vtkReebLabel* lb;

  vtkIdType A0 = n->ArcDownId;
  vtkIdType A1 = n->ArcUpId;

  vtkReebArc* a0 = this->GetArc(A0);
  vtkReebArc* a1 = this->GetArc(A1);

  a0->NodeId1 = a1->NodeId1;
  a0->ArcUpId1 = a1->ArcUpId1;

  if (a1->ArcUpId1)
    this->GetArc(a1->ArcUpId1)->ArcDwId1 = A0;

  a0->ArcDwId1 = a1->ArcDwId1;

  if (a1->ArcDwId1)
    this->GetArc(a1->ArcDwId1)->ArcUpId1 = A0;

  if (this->GetNode(a1->NodeId1)->ArcDownId == A1)
    this->GetNode(a1->NodeId1)->ArcDownId = A0;

  for (Lb = a1->LabelId0; Lb; Lb = Lnext)
  {
    lb = this->GetLabel(Lb);
    Lnext = lb->HNext;

    if (lb->VPrev)
    {
      La = lb->VPrev;
      this->GetLabel(La)->VNext = lb->VNext;
    }

    if (lb->VNext)
      this->GetLabel(lb->VNext)->VPrev = lb->VPrev;

    // delete the label...
    this->GetLabel(Lb)->HNext = ((int)-2);
    this->GetLabel(Lb)->ArcId = this->MainLabelTable.FreeZone;
    this->MainLabelTable.FreeZone = (Lb);
    --(this->MainLabelTable.Number);
  }

  // delete the arc from the graph...
  this->GetArc(A1)->LabelId1 = ((int)-2);
  this->GetArc(A1)->LabelId0 = this->MainArcTable.FreeZone;
  this->MainArcTable.FreeZone = (A1);
  --(this->MainArcTable.Number);

  // delete the node from the graph...
  this->GetNode(N)->ArcUpId = ((int)-2);
  this->GetNode(N)->ArcDownId = this->MainNodeTable.FreeZone;
  this->MainNodeTable.FreeZone = (N);
  --(this->MainNodeTable.Number);
}

//------------------------------------------------------------------------------
void vtkReebGraph::Implementation::DeepCopy(Implementation* srcG)
{
  MinimumScalarValue = srcG->MinimumScalarValue;
  MaximumScalarValue = srcG->MaximumScalarValue;

  inputMesh = srcG->inputMesh;
  inputScalarField = srcG->inputScalarField;

  ArcNumber = srcG->ArcNumber;
  NodeNumber = srcG->NodeNumber;
  LoopNumber = srcG->LoopNumber;

  ScalarField = srcG->ScalarField;

  VertexStream = srcG->VertexStream;

  if (srcG->MainArcTable.Buffer)
  {
    MainArcTable.Size = srcG->MainArcTable.Size;
    MainArcTable.Number = srcG->MainArcTable.Number;
    free(this->MainArcTable.Buffer);
    this->MainArcTable.Buffer = (vtkReebArc*)malloc(sizeof(vtkReebArc) * srcG->MainArcTable.Size);

    memcpy(this->MainArcTable.Buffer, srcG->MainArcTable.Buffer,
      sizeof(vtkReebArc) * srcG->MainArcTable.Size);
  }

  if (srcG->MainNodeTable.Buffer)
  {
    MainNodeTable.Size = srcG->MainNodeTable.Size;
    MainNodeTable.Number = srcG->MainNodeTable.Number;
    // free existing buffer
    free(this->MainNodeTable.Buffer);
    this->MainNodeTable.Buffer =
      (vtkReebNode*)malloc(sizeof(vtkReebNode) * srcG->MainNodeTable.Size);

    memcpy(this->MainNodeTable.Buffer, srcG->MainNodeTable.Buffer,
      sizeof(vtkReebNode) * srcG->MainNodeTable.Size);
  }

  if (srcG->MainLabelTable.Buffer)
  {
    MainLabelTable.Size = srcG->MainLabelTable.Size;
    MainLabelTable.Number = srcG->MainLabelTable.Number;
    free(this->MainLabelTable.Buffer);

    this->MainLabelTable.Buffer =
      (vtkReebLabel*)malloc(sizeof(vtkReebLabel) * srcG->MainLabelTable.Size);
    memcpy(this->MainLabelTable.Buffer, srcG->MainLabelTable.Buffer,
      sizeof(vtkReebLabel) * srcG->MainLabelTable.Size);
  }

  if (srcG->ArcLoopTable)
  {
    this->ArcLoopTable = (vtkIdType*)malloc(sizeof(vtkIdType) * srcG->LoopNumber);
    memcpy(this->ArcLoopTable, srcG->ArcLoopTable, sizeof(vtkIdType) * srcG->LoopNumber);
  }

  if (srcG->VertexMapSize)
  {
    this->VertexMapSize = srcG->VertexMapSize;
    this->VertexMapAllocatedSize = srcG->VertexMapAllocatedSize;
    this->VertexMap = (vtkIdType*)malloc(sizeof(vtkIdType) * this->VertexMapAllocatedSize);
    memcpy(this->VertexMap, srcG->VertexMap, sizeof(vtkIdType) * srcG->VertexMapAllocatedSize);
  }

  if (srcG->TriangleVertexMapSize)
  {
    this->TriangleVertexMapSize = srcG->TriangleVertexMapSize;
    this->TriangleVertexMapAllocatedSize = srcG->TriangleVertexMapAllocatedSize;
    this->TriangleVertexMap = (int*)malloc(sizeof(int) * this->TriangleVertexMapAllocatedSize);
    memcpy(this->TriangleVertexMap, srcG->TriangleVertexMap,
      sizeof(int) * srcG->TriangleVertexMapAllocatedSize);
  }
}

//------------------------------------------------------------------------------
vtkStandardNewMacro(vtkReebGraph);

//------------------------------------------------------------------------------
void vtkReebGraph::Implementation::SetLabel(vtkIdType arcId, vtkReebLabelTag Label)
{
  inputMesh = nullptr;

  ResizeMainLabelTable(1);

  // create a new label in the graph
  vtkIdType L = this->MainLabelTable.FreeZone;
  this->MainLabelTable.FreeZone = this->GetLabel(L)->ArcId;
  ++(this->MainLabelTable.Number);
  memset(this->GetLabel(L), 0, sizeof(vtkReebLabel));
  vtkReebLabel* l = this->GetLabel(L);

  l->HPrev = 0;
  l->HNext = 0;
  this->GetArc(arcId)->LabelId0 = L;
  this->GetArc(arcId)->LabelId1 = L;

  l->ArcId = arcId;
  l->label = Label;

  vtkIdType Lp = FindDwLabel(this->GetArc(arcId)->NodeId0, Label);
  vtkIdType Ln = FindUpLabel(this->GetArc(arcId)->NodeId1, Label);

  l->VPrev = Lp;
  if (Lp)
    this->GetLabel(Lp)->VNext = L;
  l->VNext = Ln;
  if (Ln)
    this->GetLabel(Ln)->VPrev = L;
}

//------------------------------------------------------------------------------
void vtkReebGraph::Implementation::FastArcSimplify(
  vtkIdType arcId, int vtkNotUsed(ArcNumber), vtkIdType* vtkNotUsed(arcTable))
{
  // Remove the arc which opens the loop
  vtkIdType nodeId0 = this->GetArc(arcId)->NodeId0;
  vtkIdType nodeId1 = this->GetArc(arcId)->NodeId1;

  vtkReebArc* A = this->GetArc(arcId);
  vtkReebArc* B = nullptr;
  int down, middle, up;

  if (historyOn)
  {
    if (A->ArcDwId0)
    {
      B = this->GetArc(A->ArcDwId0);

      down = this->GetNode(B->NodeId0)->VertexId;
      middle = this->GetNode(A->NodeId0)->VertexId;
      up = this->GetNode(B->NodeId1)->VertexId;

      vtkReebCancellation c;
      c.removedArcs.emplace_back(middle, up);
      c.insertedArcs.emplace_back(down, up);
      this->cancellationHistory.push_back(c);
    }
    if (A->ArcDwId1)
    {
      B = this->GetArc(A->ArcDwId1);

      down = this->GetNode(B->NodeId0)->VertexId;
      middle = this->GetNode(A->NodeId0)->VertexId;
      up = this->GetNode(A->NodeId1)->VertexId;

      vtkReebCancellation c;
      c.removedArcs.emplace_back(middle, up);
      c.insertedArcs.emplace_back(down, up);
      this->cancellationHistory.push_back(c);
    }
    if (A->ArcUpId0)
    {
      B = this->GetArc(A->ArcUpId0);

      down = this->GetNode(A->NodeId0)->VertexId;
      middle = this->GetNode(A->NodeId1)->VertexId;
      up = this->GetNode(B->NodeId1)->VertexId;

      vtkReebCancellation c;
      c.removedArcs.emplace_back(down, middle);
      c.insertedArcs.emplace_back(down, up);
      this->cancellationHistory.push_back(c);
    }
    if (A->ArcUpId1)
    {
      B = this->GetArc(A->ArcUpId1);

      down = this->GetNode(B->NodeId0)->VertexId;
      middle = this->GetNode(A->NodeId1)->VertexId;
      up = this->GetNode(B->NodeId1)->VertexId;

      vtkReebCancellation c;
      c.removedArcs.emplace_back(down, middle);
      c.insertedArcs.emplace_back(down, up);
      this->cancellationHistory.push_back(c);
    }
  }

  vtkReebGraphRemoveUpArc(this, nodeId0, arcId);
  vtkReebGraphRemoveDownArc(this, nodeId1, arcId);

  // delete the arc from the graph...
  this->GetArc(arcId)->LabelId1 = ((int)-2);
  this->GetArc(arcId)->LabelId0 = this->MainArcTable.FreeZone;
  this->MainArcTable.FreeZone = (arcId);
  --(this->MainArcTable.Number);
}

//------------------------------------------------------------------------------
vtkIdType vtkReebGraph::Implementation::FindGreater(
  vtkIdType nodeId, vtkIdType startingNodeId, vtkReebLabelTag label)
{
  if (!this->GetNode(nodeId)->IsFinalized)
    return 0;

  // base case
  if (vtkReebGraphIsHigherThan(
        this, nodeId, startingNodeId, this->GetNode(nodeId), this->GetNode(startingNodeId)))
    return nodeId;

  // iterative case
  for (vtkIdType A = this->GetNode(nodeId)->ArcUpId; A; A = this->GetArc(A)->ArcDwId0)
  {
    vtkReebArc* a = this->GetArc(A);
    vtkIdType M = this->GetArc(A)->NodeId1;
    vtkReebNode* m = this->GetNode(M);

    if (a->LabelId0 || !m->IsFinalized) // other labels or not final node
    {
      continue;
    }

    if ((M = FindGreater(M, startingNodeId, label)))
    {
      if (label)
      {
        SetLabel(A, label);
      }
      return M;
    }
  }

  return 0;
}

//------------------------------------------------------------------------------
vtkIdType vtkReebGraph::Implementation::FindLess(
  vtkIdType nodeId, vtkIdType startingNodeId, vtkReebLabelTag label)
{
  if (!this->GetNode(nodeId)->IsFinalized)
    return 0;

  // base case
  if (vtkReebGraphIsSmaller(
        this, nodeId, startingNodeId, this->GetNode(nodeId), this->GetNode(startingNodeId)))
    return nodeId;

  // iterative case
  for (vtkIdType A = this->GetNode(nodeId)->ArcDownId; A; A = this->GetArc(A)->ArcDwId1)
  {
    vtkReebArc* a = this->GetArc(A);
    vtkIdType M = this->GetArc(A)->NodeId0;
    vtkReebNode* m = this->GetNode(M);

    if (a->LabelId0 || !m->IsFinalized) // other labels or not final node
      continue;

    if ((M = FindLess(M, startingNodeId, label)))
    {
      if (label)
      {
        SetLabel(A, label);
      }
      return M;
    }
  }

  return 0;
}

//------------------------------------------------------------------------------
vtkIdType vtkReebGraph::Implementation::FindJoinNode(
  vtkIdType arcId, vtkReebLabelTag label, bool onePathOnly)
{
  vtkIdType N = this->GetArc(arcId)->NodeId1;
  vtkIdType Ret;

  if (this->GetArc(arcId)->LabelId0 || !this->GetNode(N)->IsFinalized)
    // other labels or not final node
    return 0;

  if (onePathOnly && (this->GetArc(arcId)->ArcDwId0 || this->GetArc(arcId)->ArcUpId0))
    return 0;

  // base case
  if (this->GetArc(arcId)->ArcDwId1 || this->GetArc(arcId)->ArcUpId1)
  {
    if (label)
      SetLabel(arcId, label);
    return N;
  }

  for (vtkIdType C = this->GetNode(N)->ArcUpId; C; C = this->GetArc(C)->ArcDwId0)
  {
    Ret = FindJoinNode(C, label, onePathOnly);

    if (Ret)
    {
      if (label)
        SetLabel(arcId, label);
      return Ret;
    }
  }

  return 0;
}

//------------------------------------------------------------------------------
vtkIdType vtkReebGraph::Implementation::FindSplitNode(
  vtkIdType arcId, vtkReebLabelTag label, bool onePathOnly)
{
  vtkIdType N = this->GetArc(arcId)->NodeId0;
  vtkIdType Ret;

  if (this->GetArc(arcId)->LabelId0 || !this->GetNode(N)->IsFinalized)
    // other labels or not final node
    return 0;

  if (onePathOnly && (this->GetArc(arcId)->ArcDwId1 || this->GetArc(arcId)->ArcUpId1))
    return 0;

  // base case
  if (this->GetArc(arcId)->ArcDwId0 || this->GetArc(arcId)->ArcUpId0)
  {
    if (label)
      SetLabel(arcId, label);
    return N;
  }

  // iterative case
  for (vtkIdType C = this->GetNode(N)->ArcDownId; C; C = this->GetArc(C)->ArcDwId1)
  {
    Ret = FindSplitNode(C, label, onePathOnly);

    if (Ret)
    {
      if (label)
        SetLabel(arcId, label);
      return Ret;
    }
  }

  return 0;
}

//------------------------------------------------------------------------------
vtkReebGraph::Implementation::vtkReebPath vtkReebGraph::Implementation::FindPath(
  vtkIdType arcId, double simplificationThreshold, vtkReebGraphSimplificationMetric* metric)
{
  vtkReebPath entry;
  std::priority_queue<vtkReebPath> pq;
  int size;

  vtkIdType N0 = this->GetArc(arcId)->NodeId0;
  vtkIdType N1 = this->GetArc(arcId)->NodeId1;

  char* Ntouch = nullptr;
  char* Atouch = nullptr;

  //  double simplificationValue = 0;
  if ((!inputMesh) || (!metric))
  {
    double f0 = this->GetNode(N0)->Value;
    double f1 = this->GetNode(N1)->Value;
    entry.SimplificationValue = (f1 - f0) / (this->MaximumScalarValue - this->MinimumScalarValue);
  }
  else
  {
    entry.SimplificationValue = ComputeCustomMetric(metric, this->GetArc(arcId));
  }

  // the arc itself has a good persistence
  if (simplificationThreshold && entry.SimplificationValue >= simplificationThreshold)
  {
  NOT_FOUND:
    if (Ntouch)
      free(Ntouch);
    if (Atouch)
      free(Atouch);
    vtkReebPath fake;
    memset(&fake, 0, sizeof(vtkReebPath));
    fake.SimplificationValue = 1;
    return fake;
  }

  Atouch = (char*)malloc(sizeof(char) * this->MainArcTable.Size);
  Ntouch = (char*)malloc(sizeof(char) * this->MainNodeTable.Size);
  memset(Atouch, 0, sizeof(char) * this->MainArcTable.Size);
  memset(Ntouch, 0, sizeof(char) * this->MainNodeTable.Size);

  Ntouch[N0] = 1;

  // I don't want to use the arc given by the user
  Atouch[arcId] = 1;

  entry.NodeNumber = 1;
  entry.NodeTable = new vtkIdType[1];
  entry.NodeTable[0] = N0;
  entry.ArcNumber = 0;
  entry.ArcTable = nullptr;
  pq.push(entry);

  while ((size = static_cast<int>(pq.size())))
  {
    entry = pq.top();
    pq.pop();

    int N = entry.NodeTable[entry.NodeNumber - 1];

    for (int dir = 0; dir <= 1; dir++)
    {
      for (vtkIdType A = (!dir) ? this->GetNode(N)->ArcDownId : this->GetNode(N)->ArcUpId; A;
           A = (!dir) ? this->GetArc(A)->ArcDwId1 : this->GetArc(A)->ArcDwId0)
      {
        int M = (!dir) ? (this->GetArc(A)->NodeId0) : (this->GetArc(A)->NodeId1);

        if (Atouch[A])
          continue;
        Atouch[A] = 1;

        // already used (==there is a better path to reach M)
        if (Ntouch[M])
          continue;
        Ntouch[M] = 1;

        // found!!!
        if (M == N1)
        {
          // clear all the items in the priority queue
          while (!pq.empty())
          {
            vtkReebPath aux = pq.top();
            pq.pop();
            delete aux.ArcTable;
            delete aux.NodeTable;
          }

          free(Ntouch);
          free(Atouch);

          vtkIdType* tmp = new vtkIdType[entry.NodeNumber + 1];
          memcpy(tmp, entry.NodeTable, sizeof(vtkIdType) * entry.NodeNumber);
          tmp[entry.NodeNumber] = N1;
          delete[] entry.NodeTable;
          entry.NodeTable = tmp;
          entry.NodeNumber++;
          return entry;
        }

        if ((!inputMesh) || (!metric))
        {
          entry.SimplificationValue += vtkReebGraphGetArcPersistence(this, this->GetArc(A));
        }
        else
        {
          entry.SimplificationValue += ComputeCustomMetric(metric, this->GetArc(A));
        }

        // The loop persistence is greater than functionScale
        if (simplificationThreshold && entry.SimplificationValue >= simplificationThreshold)
          continue;

        vtkReebPath newentry;
        newentry.SimplificationValue = entry.SimplificationValue;
        newentry.ArcNumber = entry.ArcNumber + 1;
        newentry.ArcTable = new vtkIdType[newentry.ArcNumber];
        newentry.NodeNumber = entry.NodeNumber + 1;
        newentry.NodeTable = new vtkIdType[newentry.NodeNumber];
        if (entry.ArcNumber)
          memcpy(newentry.ArcTable, entry.ArcTable, sizeof(vtkIdType) * entry.ArcNumber);
        if (entry.NodeNumber)
          memcpy(newentry.NodeTable, entry.NodeTable, sizeof(vtkIdType) * entry.NodeNumber);

        newentry.ArcTable[entry.ArcNumber] = A;
        newentry.NodeTable[entry.NodeNumber] = M;
        pq.push(newentry);
      }
    }

    // finished with this entry
    delete entry.ArcTable;
    delete[] entry.NodeTable;
  }

  goto NOT_FOUND;
}

//------------------------------------------------------------------------------
int vtkReebGraph::Implementation::SimplifyLoops(
  double simplificationThreshold, vtkReebGraphSimplificationMetric* simplificationMetric)
{

  if (!simplificationThreshold)
    return 0;

  // refresh information about ArcLoopTable
  this->FindLoops();

  int NumSimplified = 0;

  for (int n = 0; n < this->LoopNumber; n++)
  {
    int A = this->ArcLoopTable[n];

    if (this->GetArc(A)->LabelId1 == ((int)-2))
      continue;

    double simplificationValue = 0;
    if ((!inputMesh) || (!simplificationMetric))
    {
      vtkIdType N0 = this->GetArc(A)->NodeId0;
      vtkIdType N1 = this->GetArc(A)->NodeId1;
      double f0 = this->GetNode(N0)->Value;
      double f1 = this->GetNode(N1)->Value;
      simplificationValue = (f1 - f0) / (this->MaximumScalarValue - this->MinimumScalarValue);
    }
    else
    {
      simplificationValue = ComputeCustomMetric(simplificationMetric, this->GetArc(A));
    }

    if (simplificationValue >= simplificationThreshold)
      continue;

    vtkReebPath entry =
      this->FindPath(this->ArcLoopTable[n], simplificationThreshold, simplificationMetric);

    // too high for persistence
    if (entry.SimplificationValue >= simplificationThreshold)
      continue;

    // distribute its bucket to the loop and delete the arc
    this->FastArcSimplify(ArcLoopTable[n], entry.ArcNumber, entry.ArcTable);
    delete entry.ArcTable;
    delete entry.NodeTable;

    ++NumSimplified;
    CommitSimplification();
  }

  // check for regular points
  for (int N = 1; N < this->MainNodeTable.Size; N++)
  {
    if (this->GetNode(N)->ArcUpId == ((int)-2))
      continue;

    if (this->GetNode(N)->ArcDownId == 0 && this->GetNode(N)->ArcUpId == 0)
    {
      // delete the node from the graph...
      this->GetNode(N)->ArcUpId = ((int)-2);
      this->GetNode(N)->ArcDownId = this->MainNodeTable.FreeZone;
      this->MainNodeTable.FreeZone = (N);
      --(this->MainNodeTable.Number);
    }

    else if (vtkReebGraphIsRegular(this, this->GetNode(N)))
    {
      if (historyOn)
      {
        vtkReebNode* n = this->GetNode(N);

        vtkIdType A0 = n->ArcDownId;
        vtkIdType A1 = n->ArcUpId;

        vtkReebArc* a0 = this->GetArc(A0);
        vtkReebArc* a1 = this->GetArc(A1);
        vtkReebNode* downN = this->GetNode(a0->NodeId0);
        vtkReebNode* upN = this->GetNode(a1->NodeId1);

        int down, middle, up;
        down = downN->VertexId;
        middle = n->VertexId;
        up = upN->VertexId;

        vtkReebCancellation c;
        c.removedArcs.emplace_back(down, middle);
        c.removedArcs.emplace_back(middle, up);
        c.insertedArcs.emplace_back(down, up);

        this->cancellationHistory.push_back(c);
      }
      EndVertex(N);
    }
  }

  this->RemovedLoopNumber = NumSimplified;

  return NumSimplified;
}

//------------------------------------------------------------------------------
double vtkReebGraph::Implementation::ComputeCustomMetric(
  vtkReebGraphSimplificationMetric* simplificationMetric, vtkReebArc* a)
{
  int edgeId = -1, start = -1, end = -1;

  vtkDataArray* vertexInfo =
    vtkArrayDownCast<vtkDataArray>(this->Parent->GetVertexData()->GetAbstractArray("Vertex Ids"));
  if (!vertexInfo)
    return vtkReebGraphGetArcPersistence(this, a);

  vtkVariantArray* edgeInfo =
    vtkArrayDownCast<vtkVariantArray>(this->Parent->GetEdgeData()->GetAbstractArray("Vertex Ids"));
  if (!edgeInfo)
    return vtkReebGraphGetArcPersistence(this, a);

  vtkEdgeListIterator* eIt = vtkEdgeListIterator::New();
  this->Parent->GetEdges(eIt);

  do
  {
    vtkEdgeType e = eIt->Next();
    if (((*(vertexInfo->GetTuple(e.Source))) == GetNodeVertexId(a->NodeId0)) &&
      ((*(vertexInfo->GetTuple(e.Target))) == GetNodeVertexId(a->NodeId1)))
    {
      edgeId = e.Id;
      start = static_cast<int>(*(vertexInfo->GetTuple(e.Source)));
      end = static_cast<int>(*(vertexInfo->GetTuple(e.Target)));
      break;
    }
  } while (eIt->HasNext());
  eIt->Delete();

  vtkAbstractArray* vertexList = edgeInfo->GetPointer(edgeId)->ToArray();

  return simplificationMetric->ComputeMetric(inputMesh, inputScalarField, start, vertexList, end);
}

//------------------------------------------------------------------------------
int vtkReebGraph::Implementation::SimplifyBranches(
  double simplificationThreshold, vtkReebGraphSimplificationMetric* simplificationMetric)
{
  static const vtkReebLabelTag RouteOld = 100;
  static const vtkReebLabelTag RouteNew = 200;
  int nstack, mstack = 0;
  int* stack = nullptr;

  if (!simplificationThreshold)
    return 0;

  int nsimp = 0;
  int cont = 0;
  const int step = 10000;
  bool redo;

  vtkDataSet* input = inputMesh;

REDO:
  nstack = 0;
  redo = false;

  for (int N = 1; N < this->MainNodeTable.Size; ++N)
  {
    if (this->GetNode(N)->ArcUpId == ((int)-2))
      continue;

    vtkReebNode* n = this->GetNode(N);

    // simplify atomic nodes
    if (!n->ArcDownId && !n->ArcUpId)
    {
      // delete the node from the graph...
      this->GetNode(N)->ArcUpId = ((int)-2);
      this->GetNode(N)->ArcDownId = this->MainNodeTable.FreeZone;
      this->MainNodeTable.FreeZone = (N);
      --(this->MainNodeTable.Number);
    }
    else if (!n->ArcDownId)
    {
      // insert into stack branches to simplify
      for (int A_ = n->ArcUpId; A_; A_ = this->GetArc(A_)->ArcDwId0)
      {
        vtkReebArc* a_ = this->GetArc(A_);
        if ((!inputMesh) || (!simplificationMetric))
        {
          if (vtkReebGraphGetArcPersistence(this, a_) < simplificationThreshold)
          {
            vtkReebGraphStackPush(A_);
          }
        }
        else
        {
          if (this->ComputeCustomMetric(simplificationMetric, a_) < simplificationThreshold)
          {
            vtkReebGraphStackPush(A_);
          }
        }
      }
    }
    else if (!n->ArcUpId)
    {
      // insert into stack branches to simplify
      for (vtkIdType A_ = n->ArcDownId; A_; A_ = this->GetArc(A_)->ArcDwId1)
      {
        vtkReebArc* a_ = this->GetArc(A_);
        if (vtkReebGraphGetArcPersistence(this, a_) < simplificationThreshold)
        {
          vtkReebGraphStackPush(A_);
        }
      }
    }
  }

  while (vtkReebGraphStackSize())
  {
    int A = vtkReebGraphStackTop();
    vtkReebGraphStackPop();

    if (!--cont)
    {
      cont = step;
    }

    if (this->GetArc(A)->LabelId1 == ((int)-2))
      continue;

    cont++;

    vtkReebArc* arc = this->GetArc(A);

    int N = arc->NodeId0;
    int M = arc->NodeId1;

    if (this->GetNode(N)->ArcDownId && this->GetNode(M)->ArcUpId)
      continue;

    double persistence = vtkReebGraphGetArcPersistence(this, arc);

    // is the actual persistence (in percentage) greater than the applied filter?
    if (persistence >= simplificationThreshold)
      continue;

    int Mdown = 0, Nup = 0, Ndown = 0, Mup = 0;

    // get the 'down' degree for M
    for (vtkIdType A_ = this->GetNode(M)->ArcDownId; A_; A_ = this->GetArc(A_)->ArcDwId1)
      ++Mdown;

    // Get the 'up' degree for N
    for (vtkIdType A_ = this->GetNode(N)->ArcUpId; A_; A_ = this->GetArc(A_)->ArcDwId0)
      ++Nup;

    // get the 'down' degree for N
    for (vtkIdType A_ = this->GetNode(N)->ArcDownId; A_; A_ = this->GetArc(A_)->ArcDwId1)
      ++Ndown;

    // get the 'up' degree for M
    for (vtkIdType A_ = this->GetNode(M)->ArcUpId; A_; A_ = this->GetArc(A_)->ArcDwId0)
      ++Mup;

    // isolated arc
    if (!Ndown && Nup == 1 && Mdown == 1 && !Mup)
    {
      vtkReebGraphRemoveUpArc(this, N, A);
      vtkReebGraphRemoveDownArc(this, M, A);

      // delete the arc from the graph...
      this->GetArc(A)->LabelId1 = ((int)-2);
      this->GetArc(A)->LabelId0 = this->MainArcTable.FreeZone;
      this->MainArcTable.FreeZone = (A);
      --(this->MainArcTable.Number);

      if (!(this->GetNode(N)->ArcUpId == ((int)-2)) &&
        vtkReebGraphIsRegular(this, this->GetNode(N)))
      {
        EndVertex(N);
      }
      if (!(this->GetNode(M)->ArcUpId == ((int)-2)) &&
        vtkReebGraphIsRegular(this, this->GetNode(M)))
      {
        EndVertex(M);
      }

      nsimp++;
      redo = true;
      continue;
    }

    int Down = 0;
    int Up = 0;

    bool simplified = false;

    // M is a maximum
    if (!Mup)
    {
      if ((Down = FindSplitNode(A, RouteOld)))
      {
        if ((Up = FindGreater(Down, M, RouteNew)))
        {
          SetLabel(AddArc(M, Up), RouteOld);
          Collapse(Down, Up, RouteOld, RouteNew);
          simplified = true;
        }
        else
        {
          this->SimplifyLabels(Down);
        }
      }
    }

    // N is a minimum
    if (!simplified && !Ndown)
    {
      if ((Up = FindJoinNode(A, RouteOld)))
      {
        if ((Down = FindLess(Up, N, RouteNew)))
        {
          SetLabel(AddArc(Down, N), RouteOld);
          Collapse(Down, Up, RouteOld, RouteNew);
          simplified = true;
        }
        else
        {
          this->SimplifyLabels(Up);
        }
      }
    }

    if (simplified)
    {
      if (!(this->GetNode(Down)->ArcUpId == ((int)-2)))
      {
        this->SimplifyLabels(Down);

        if (!this->GetNode(Down)->ArcDownId) // minimum
        {
          for (vtkIdType A_ = this->GetNode(Down)->ArcUpId; A_; A_ = this->GetArc(A_)->ArcDwId0)
          {
            vtkReebArc* a_ = this->GetArc(A_);
            if ((!inputMesh) || (!simplificationMetric))
            {
              if (vtkReebGraphGetArcPersistence(this, a_) < simplificationThreshold)
              {
                vtkReebGraphStackPush(A_);
              }
            }
            else
            {
              if (this->ComputeCustomMetric(simplificationMetric, a_) < simplificationThreshold)
              {
                vtkReebGraphStackPush(A_);
              }
            }
          }
        }
      }

      if (!(this->GetNode(Up)->ArcUpId == ((int)-2)))
      {
        this->SimplifyLabels(Up);

        if (!this->GetNode(Up)->ArcUpId)
        {
          for (vtkIdType A_ = this->GetNode(Up)->ArcDownId; A_; A_ = this->GetArc(A_)->ArcDwId1)
          {
            vtkReebArc* a_ = this->GetArc(A_);
            if ((!inputMesh) || (!simplificationMetric))
            {
              if (vtkReebGraphGetArcPersistence(this, a_) < simplificationThreshold)
              {
                vtkReebGraphStackPush(A_);
              }
            }
            else
            {
              if (this->ComputeCustomMetric(simplificationMetric, a_) < simplificationThreshold)
              {
                vtkReebGraphStackPush(A_);
              }
            }
          }
        }
      }

      nsimp++;
      redo = true;
    }
    CommitSimplification();
  } // while

  if (redo)
    goto REDO;

  free(stack);

  inputMesh = input;

  return nsimp;
}

//------------------------------------------------------------------------------
void vtkReebGraph::Implementation::ResizeMainNodeTable(int newSize)
{
  int oldsize, i;

  if ((this->MainNodeTable.Size - this->MainNodeTable.Number) < newSize)
  {
    oldsize = this->MainNodeTable.Size;

    if (!this->MainNodeTable.Size)
      this->MainNodeTable.Size = newSize;
    while ((this->MainNodeTable.Size - this->MainNodeTable.Number) < newSize)
      this->MainNodeTable.Size <<= 1;

    this->MainNodeTable.Buffer = (vtkReebNode*)realloc(
      this->MainNodeTable.Buffer, sizeof(vtkReebNode) * this->MainNodeTable.Size);

    for (i = oldsize; i < this->MainNodeTable.Size - 1; i++)
    {
      this->GetNode(i)->ArcDownId = i + 1;
      this->GetNode(i)->ArcUpId = ((int)-2);
    }

    this->GetNode(i)->ArcDownId = this->MainNodeTable.FreeZone;
    this->GetNode(i)->ArcUpId = ((int)-2);
    this->MainNodeTable.FreeZone = oldsize;
  }
}

//------------------------------------------------------------------------------
int vtkReebGraph::Implementation::CommitSimplification()
{
  // now re-construct the graph with projected deg-2 nodes.
  std::vector<std::pair<std::pair<int, int>, std::vector<int>>> before, after;

  vtkEdgeListIterator* eIt = vtkEdgeListIterator::New();
  this->Parent->GetEdges(eIt);
  vtkVariantArray* edgeInfo =
    vtkArrayDownCast<vtkVariantArray>(this->Parent->GetEdgeData()->GetAbstractArray("Vertex Ids"));
  vtkDataArray* vertexInfo = this->Parent->GetVertexData()->GetArray("Vertex Ids");

  // avoids double projection
  int vertexNumber = vertexInfo->GetNumberOfTuples();
  std::vector<bool> segmentedVertices;

  do
  {
    std::pair<std::pair<int, int>, std::vector<int>> superArc;

    vtkEdgeType e = eIt->Next();
    vtkAbstractArray* vertexList = edgeInfo->GetPointer(e.Id)->ToArray();

    superArc.first.first = (int)*(vertexInfo->GetTuple(e.Source));
    superArc.first.second = (int)*(vertexInfo->GetTuple(e.Target));

    superArc.second.resize(vertexList->GetNumberOfTuples());
    vertexNumber += vertexList->GetNumberOfTuples();
    for (unsigned int i = 0; i < superArc.second.size(); i++)
      superArc.second[i] = vertexList->GetVariantValue(i).ToInt();

    before.push_back(superArc);
  } while (eIt->HasNext());

  segmentedVertices.resize(vertexNumber);
  for (unsigned int i = 0; i < segmentedVertices.size(); i++)
    segmentedVertices[i] = false;

  vtkIdType prevArcId = -1, arcId = 0;
  while (arcId != prevArcId)
  {
    prevArcId = arcId;
    arcId = GetPreviousArcId();
  }
  prevArcId = -1;

  while (prevArcId != arcId)
  {
    if (this->GetArc(arcId))
    {
      int down, up;
      down = this->GetNode((this->GetArc(arcId))->NodeId0)->VertexId;
      up = this->GetNode((this->GetArc(arcId))->NodeId1)->VertexId;

      std::pair<std::pair<int, int>, std::vector<int>> superArc;

      superArc.first.first = down;
      superArc.first.second = up;

      after.push_back(superArc);
    }
    prevArcId = arcId;
    arcId = GetNextArcId();
  }

  std::pair<int, int> destinationArc;
  std::map<int, bool> processedOutputArcs;

  // now map the unsimplified arcs onto the simplified ones
  for (unsigned int i = 0; i < before.size(); i++)
  {
    std::vector<int> simplifiedCriticalNodes;
    destinationArc = before[i].first;
    for (unsigned int j = 0; j < this->cancellationHistory.size(); j++)
    {
      for (unsigned int k = 0; k < this->cancellationHistory[j].removedArcs.size(); k++)
      {
        if ((destinationArc.first == this->cancellationHistory[j].removedArcs[k].first) &&
          (destinationArc.second == this->cancellationHistory[j].removedArcs[k].second))
        {
          // the arc has been involved in a cancellation
          destinationArc = this->cancellationHistory[j].insertedArcs[0];

          if (this->cancellationHistory[j].removedArcs.size() > 1)
          {
            if (((this->cancellationHistory[j].removedArcs[0].first == destinationArc.first) &&
                  (this->cancellationHistory[j].removedArcs[1].second == destinationArc.second)) ||
              ((this->cancellationHistory[j].removedArcs[1].first == destinationArc.first) &&
                (this->cancellationHistory[j].removedArcs[0].second == destinationArc.second)))
            {
              for (unsigned int l = 0; l < this->cancellationHistory[j].removedArcs.size(); l++)
              {
                if ((this->cancellationHistory[j].removedArcs[l].first != destinationArc.first) &&
                  (this->cancellationHistory[j].removedArcs[l].first != destinationArc.second))
                {
                  // this critical node will become a degree two node, let's
                  // remember it
                  simplifiedCriticalNodes.push_back(
                    this->cancellationHistory[j].removedArcs[l].first);
                }
                if ((this->cancellationHistory[j].removedArcs[l].second != destinationArc.first) &&
                  (this->cancellationHistory[j].removedArcs[l].second != destinationArc.second))
                {
                  // same thing as above
                  simplifiedCriticalNodes.push_back(
                    this->cancellationHistory[j].removedArcs[l].second);
                }
              }
            }
          }
        }
      }
    }

    // at this point the deg2-nodes are in before[i].second

    // now find the projection in the simplified graph
    for (unsigned int j = 0; j < after.size(); j++)
    {
      if (destinationArc == after[j].first)
      {
        std::map<int, bool>::iterator aIt;
        aIt = processedOutputArcs.find(j);
        if (aIt == processedOutputArcs.end())
        {
          if (before[i].first == destinationArc)
          {
            // non-simplified arc
            processedOutputArcs[j] = true;
            after[j].second = before[i].second;
          }
          if (before[i].first != destinationArc)
          {
            // adding content of before[i].second to after[j].second
            for (unsigned int k = 0; k < before[i].second.size(); k++)
            {
              if (!segmentedVertices[before[i].second[k]])
              {
                after[j].second.push_back(before[i].second[k]);
                segmentedVertices[before[i].second[k]] = true;
              }
            }
          }
          for (unsigned int k = 0; k < simplifiedCriticalNodes.size(); k++)
          {
            if (!segmentedVertices[simplifiedCriticalNodes[k]])
            {
              after[j].second.push_back(simplifiedCriticalNodes[k]);
              segmentedVertices[simplifiedCriticalNodes[k]] = true;
            }
          }
          break;
        }
      }
    }
  }
  // ensure the sorting on the arcs
  for (unsigned int i = 0; i < after.size(); i++)
  {
    std::vector<std::pair<int, double>> scalarValues;
    for (unsigned int j = 0; j < after[i].second.size(); j++)
    {
      std::pair<int, double> scalarVertex;
      scalarVertex.first = after[i].second[j];
      std::map<int, double>::iterator sIt;
      sIt = ScalarField.find(scalarVertex.first);
      if (sIt != ScalarField.end())
      {
        scalarVertex.second = sIt->second;
        scalarValues.push_back(scalarVertex);
      }
    }
    std::sort(scalarValues.begin(), scalarValues.end(), vtkReebGraphVertexSoS);
    for (unsigned int j = 0; j < after[i].second.size(); j++)
      after[i].second[j] = scalarValues[j].first;
  }

  // now construct the vtkMutableDirectedGraph
  // first, clean up the current graph
  while (this->Parent->GetNumberOfEdges())
    this->Parent->RemoveEdge(0);
  while (this->Parent->GetNumberOfVertices())
    this->Parent->RemoveVertex(0);

  this->Parent->GetVertexData()->RemoveArray("Vertex Ids");
  this->Parent->GetEdgeData()->RemoveArray("Vertex Ids");

  vtkIdType prevNodeId = -1, nodeId = 0;
  while (prevNodeId != nodeId)
  {
    prevNodeId = nodeId;
    nodeId = GetPreviousNodeId();
  }
  prevNodeId = -1;

  vtkVariantArray* vertexProperties = vtkVariantArray::New();
  vertexProperties->SetNumberOfValues(1);

  vtkIdTypeArray* vertexIds = vtkIdTypeArray::New();
  vertexIds->SetName("Vertex Ids");
  this->Parent->GetVertexData()->AddArray(vertexIds);

  std::map<int, int> vMap;
  int vIt = 0;
  while (prevNodeId != nodeId)
  {
    vtkIdType nodeVertexId = GetNodeVertexId(nodeId);
    vMap[nodeVertexId] = vIt;
    vertexProperties->SetValue(0, nodeVertexId);
    this->Parent->AddVertex(vertexProperties);

    prevNodeId = nodeId;
    nodeId = GetNextNodeId();
    vIt++;
  }
  vertexIds->Delete();
  vertexProperties->Delete();

  vtkVariantArray* deg2NodeIds = vtkVariantArray::New();
  deg2NodeIds->SetName("Vertex Ids");
  this->Parent->GetEdgeData()->AddArray(deg2NodeIds);

  for (unsigned int i = 0; i < after.size(); i++)
  {
    std::map<int, int>::iterator downIt, upIt;
    downIt = vMap.find(after[i].first.first);
    upIt = vMap.find(after[i].first.second);

    if ((downIt != vMap.end()) && (upIt != vMap.end()))
    {
      vtkVariantArray* edgeProperties = vtkVariantArray::New();
      vtkIdTypeArray* vertexList = vtkIdTypeArray::New();
      vertexList->SetNumberOfValues(static_cast<vtkIdType>(after[i].second.size()));
      for (unsigned int j = 0; j < after[i].second.size(); j++)
        vertexList->SetValue(j, after[i].second[j]);
      edgeProperties->SetNumberOfValues(1);
      edgeProperties->SetValue(0, vertexList);
      this->Parent->AddEdge(downIt->second, upIt->second, edgeProperties);
      vertexList->Delete();
      edgeProperties->Delete();
    }
  }
  deg2NodeIds->Delete();

  this->cancellationHistory.clear();

  return 0;
}

//------------------------------------------------------------------------------
int vtkReebGraph::Simplify(
  double simplificationThreshold, vtkReebGraphSimplificationMetric* simplificationMetric)
{
  int deletionNumber = 0;

  this->Storage->cancellationHistory.clear();
  this->Storage->historyOn = true;

  this->Storage->ArcNumber = 0;
  this->Storage->NodeNumber = 0;

  deletionNumber = this->Storage->SimplifyBranches(simplificationThreshold, simplificationMetric) +
    this->Storage->SimplifyLoops(simplificationThreshold, simplificationMetric) +
    this->Storage->SimplifyBranches(simplificationThreshold, simplificationMetric);

  this->Storage->historyOn = false;

  return deletionNumber;
}

//------------------------------------------------------------------------------
void vtkReebGraph::Implementation::FlushLabels()
{
  for (int A = 1; A < this->MainArcTable.Size; A++)
  {
    if (!(this->GetArc(A)->LabelId1 == ((int)-2)))
      this->GetArc(A)->LabelId0 = this->GetArc(A)->LabelId1 = 0;
  }

  if (this->MainLabelTable.Buffer)
  {
    free(this->MainLabelTable.Buffer);
  }

  this->MainLabelTable.Buffer = (vtkReebLabel*)malloc(sizeof(vtkReebLabel) * 2);
  this->MainLabelTable.Size = 2;
  this->MainLabelTable.Number = 1;
  this->MainLabelTable.FreeZone = 1;
  this->GetLabel(1)->HNext = ((int)-2);
  this->GetLabel(1)->ArcId = 0;
}

//------------------------------------------------------------------------------
void vtkReebGraph::DeepCopy(vtkDataObject* src)
{

  vtkReebGraph* srcG = vtkReebGraph::SafeDownCast(src);

  if (srcG)
  {
    this->Storage->DeepCopy(srcG->Storage);
  }

  vtkMutableDirectedGraph::DeepCopy(srcG);
}

//------------------------------------------------------------------------------
void vtkReebGraph::Set(vtkMutableDirectedGraph* g)
{
  vtkMutableDirectedGraph::DeepCopy(g);
}

//------------------------------------------------------------------------------
void vtkReebGraph::CloseStream()
{

  vtkIdType prevArcId = -1, arcId = 0;
  while (arcId != prevArcId)
  {
    prevArcId = arcId;
    arcId = this->Storage->GetPreviousArcId();
  }
  prevArcId = -1;

  // loop over the arcs and build the local adjacency map

  // vertex -> (down vertices, up vertices)
  std::map<int, std::pair<std::vector<int>, std::vector<int>>> localAdjacency;
  while (prevArcId != arcId)
  {
    vtkIdType downVertexId, upVertexId;
    downVertexId = this->Storage->GetNode((this->Storage->GetArc(arcId))->NodeId0)->VertexId;
    upVertexId = this->Storage->GetNode((this->Storage->GetArc(arcId))->NodeId1)->VertexId;

    std::map<int, std::pair<std::vector<int>, std::vector<int>>>::iterator aIt;

    // lookUp for the down vertex
    aIt = localAdjacency.find(downVertexId);
    if (aIt == localAdjacency.end())
    {
      std::pair<std::vector<int>, std::vector<int>> adjacencyItem;
      adjacencyItem.second.push_back(upVertexId);
      localAdjacency[downVertexId] = adjacencyItem;
    }
    else
    {
      aIt->second.second.push_back(upVertexId);
    }

    // same thing for the up vertex
    aIt = localAdjacency.find(upVertexId);
    if (aIt == localAdjacency.end())
    {
      std::pair<std::vector<int>, std::vector<int>> adjacencyItem;
      adjacencyItem.first.push_back(downVertexId);
      localAdjacency[upVertexId] = adjacencyItem;
    }
    else
    {
      aIt->second.first.push_back(downVertexId);
    }

    prevArcId = arcId;
    arcId = this->Storage->GetNextArcId();
  }

  // now build the super-arcs with deg-2 nodes

  // <vertex,vertex>,<vertex list> (arc, deg2 node list)
  std::vector<std::pair<std::pair<int, int>, std::vector<int>>> globalAdjacency;

  std::map<int, std::pair<std::vector<int>, std::vector<int>>>::iterator aIt;
  aIt = localAdjacency.begin();
  do
  {
    if (!((aIt->second.first.size() == 1) && (aIt->second.second.size() == 1)))
    {
      // not a deg-2 node
      if (!aIt->second.second.empty())
      {
        // start the sweep up
        for (unsigned int i = 0; i < aIt->second.second.size(); i++)
        {
          std::vector<int> deg2List;
          std::map<int, std::pair<std::vector<int>, std::vector<int>>>::iterator nextIt;

          nextIt = localAdjacency.find(aIt->second.second[i]);
          while ((nextIt->second.first.size() == 1) && (nextIt->second.second.size() == 1))
          {
            deg2List.push_back(nextIt->first);
            nextIt = localAdjacency.find(nextIt->second.second[0]);
          }
          globalAdjacency.emplace_back(std::pair<int, int>(aIt->first, nextIt->first), deg2List);
        }
      }
    }
    ++aIt;
  } while (aIt != localAdjacency.end());

  // now cleanup the internal representation
  int nmyend = 0;
  for (vtkIdType N = 1; N < this->Storage->MainNodeTable.Size; N++)
  {
    // clear the node
    if (this->Storage->GetNode(N)->ArcUpId == ((int)-2))
      continue;

    vtkReebGraph::Implementation::vtkReebNode* n = this->Storage->GetNode(N);

    if (!n->IsFinalized)
    {
      nmyend++;
      this->Storage->EndVertex(N);
    }
  }

  this->Storage->FlushLabels();

  // now construct the actual graph
  vtkIdType prevNodeId = -1, nodeId = 0;
  while (prevNodeId != nodeId)
  {
    prevNodeId = nodeId;
    nodeId = this->Storage->GetPreviousNodeId();
  }
  prevNodeId = -1;

  vtkVariantArray* vertexProperties = vtkVariantArray::New();
  vertexProperties->SetNumberOfValues(1);

  vtkIdTypeArray* vertexIds = vtkIdTypeArray::New();
  vertexIds->SetName("Vertex Ids");
  GetVertexData()->AddArray(vertexIds);

  std::map<int, int> vMap;
  int vIt = 0;

  while (prevNodeId != nodeId)
  {
    vtkIdType nodeVertexId = this->Storage->GetNodeVertexId(nodeId);
    vMap[nodeVertexId] = vIt;
    vertexProperties->SetValue(0, nodeVertexId);
    AddVertex(vertexProperties);

    prevNodeId = nodeId;
    nodeId = this->Storage->GetNextNodeId();
    vIt++;
  }
  vertexIds->Delete();
  vertexProperties->Delete();

  vtkVariantArray* deg2NodeIds = vtkVariantArray::New();
  deg2NodeIds->SetName("Vertex Ids");
  GetEdgeData()->AddArray(deg2NodeIds);

  for (unsigned int i = 0; i < globalAdjacency.size(); i++)
  {
    std::map<int, int>::iterator downIt, upIt;
    downIt = vMap.find(globalAdjacency[i].first.first);
    upIt = vMap.find(globalAdjacency[i].first.second);

    if ((downIt != vMap.end()) && (upIt != vMap.end()))
    {
      vtkVariantArray* edgeProperties = vtkVariantArray::New();
      vtkIdTypeArray* vertexList = vtkIdTypeArray::New();
      vertexList->SetNumberOfValues(static_cast<vtkIdType>(globalAdjacency[i].second.size()));
      for (unsigned int j = 0; j < globalAdjacency[i].second.size(); j++)
        vertexList->SetValue(j, globalAdjacency[i].second[j]);
      edgeProperties->SetNumberOfValues(1);
      edgeProperties->SetValue(0, vertexList);
      AddEdge(downIt->second, upIt->second, edgeProperties);
      vertexList->Delete();
      edgeProperties->Delete();
    }
  }
  deg2NodeIds->Delete();
}

//------------------------------------------------------------------------------
vtkReebGraph::vtkReebGraph()
{
  this->Storage = new vtkReebGraph::Implementation;
  this->Storage->Parent = this;
}

//------------------------------------------------------------------------------
vtkReebGraph::~vtkReebGraph()
{
  delete this->Storage;
}

//------------------------------------------------------------------------------
void vtkReebGraph::PrintSelf(ostream& os, vtkIndent indent)
{
  Superclass::PrintSelf(os, indent);
  os << indent << "Reeb graph general statistics:" << endl;
  os << indent << indent << "Number Of Node(s): " << this->Storage->GetNumberOfNodes() << endl;
  os << indent << indent << "Number Of Arc(s): " << this->Storage->GetNumberOfArcs() << endl;
  os << indent << indent
     << "Number Of Connected Component(s): " << this->Storage->GetNumberOfConnectedComponents()
     << endl;
  os << indent << indent << "Number Of Loop(s): " << this->Storage->GetNumberOfLoops() << endl;
}

void vtkReebGraph::PrintNodeData(ostream& os, vtkIndent indent)
{
  vtkIdType arcId = 0, nodeId = 0;
  os << indent << "Node Data:" << endl;
  vtkIdType prevNodeId = -1;

  // roll back to the beginning of the list
  while (prevNodeId != nodeId)
  {
    prevNodeId = nodeId;
    nodeId = this->Storage->GetPreviousNodeId();
  }
  prevNodeId = -1;

  while (prevNodeId != nodeId)
  {
    prevNodeId = nodeId;
    vtkIdList* downArcIdList = vtkIdList::New();
    vtkIdList* upArcIdList = vtkIdList::New();

    this->Storage->GetNodeDownArcIds(nodeId, downArcIdList);
    this->Storage->GetNodeUpArcIds(nodeId, upArcIdList);

    cout << indent << indent << "Node " << nodeId << ":" << endl;
    cout << indent << indent << indent;
    cout << "Vert: " << this->Storage->GetNodeVertexId(nodeId);
    cout << ", Val: " << this->Storage->GetNodeScalarValue(nodeId);
    cout << ", DwA:";
    for (vtkIdType i = 0; i < downArcIdList->GetNumberOfIds(); i++)
      cout << " " << this->Storage->GetArcDownNodeId(downArcIdList->GetId(i));
    cout << ", UpA:";
    for (vtkIdType i = 0; i < upArcIdList->GetNumberOfIds(); i++)
      cout << " " << this->Storage->GetArcUpNodeId(upArcIdList->GetId(i));
    cout << endl;

    downArcIdList->Delete();
    upArcIdList->Delete();
    nodeId = this->Storage->GetNextNodeId();
  }

  os << indent << "Arc Data:" << endl;
  vtkIdType prevArcId = -1;
  arcId = 0;

  // roll back to the beginning of the list
  while (prevArcId != arcId)
  {
    prevArcId = arcId;
    arcId = this->Storage->GetPreviousArcId();
  }
  prevArcId = -1;

  while (prevArcId != arcId)
  {
    prevArcId = arcId;
    cout << indent << indent << "Arc " << arcId << ":" << endl;
    cout << indent << indent << indent;
    cout << "Down: " << this->Storage->GetArcDownNodeId(arcId);
    cout << ", Up: " << this->Storage->GetArcUpNodeId(arcId);
    cout << ", Persistence: "
         << this->Storage->GetNodeScalarValue(this->Storage->GetArcUpNodeId(arcId)) -
        this->Storage->GetNodeScalarValue(this->Storage->GetArcDownNodeId(arcId));
    cout << endl;
    arcId = this->Storage->GetNextArcId();
  }
}

//------------------------------------------------------------------------------
void vtkReebGraph::Implementation::GetNodeDownArcIds(vtkIdType nodeId, vtkIdList* arcIdList)
{
  vtkIdType i = 0;

  if (!arcIdList)
    return;

  arcIdList->Reset();

  for (vtkIdType arcId = this->GetNode(nodeId)->ArcDownId; arcId;
       arcId = this->GetArc(arcId)->ArcDwId1)
  {
    arcIdList->InsertId(i, arcId);
    i++;
  }
}

//------------------------------------------------------------------------------
void vtkReebGraph::Implementation::GetNodeUpArcIds(vtkIdType nodeId, vtkIdList* arcIdList)
{
  vtkIdType i = 0;

  if (!arcIdList)
    return;

  for (vtkIdType arcId = this->GetNode(nodeId)->ArcUpId; arcId;
       arcId = this->GetArc(arcId)->ArcDwId0)
  {
    arcIdList->InsertId(i, arcId);
    i++;
  }
}

//------------------------------------------------------------------------------
void vtkReebGraph::Implementation::FindLoops()
{

  if (this->ArcLoopTable)
  {
    free(this->ArcLoopTable);
    this->ArcLoopTable = nullptr;
    this->LoopNumber = 0;
  }

  this->ConnectedComponentNumber = 0;

  int nstack = 0, mstack = 0;
  int* stack = nullptr;

  char* Ntouch = (char*)malloc(sizeof(char) * this->MainNodeTable.Size);
  char* Atouch = (char*)malloc(sizeof(char) * this->MainArcTable.Size);

  memset(Ntouch, 0, sizeof(char) * this->MainNodeTable.Size);

  for (int Node = 1; Node < this->MainNodeTable.Size; Node++)
  {
    // check that the node is clear
    if (this->GetNode(Node)->ArcUpId == ((int)-2))
      continue;

    if (!Ntouch[Node])
    {
      ++(this->ConnectedComponentNumber);

      memset(Atouch, 0, sizeof(bool) * this->MainArcTable.Size);

      Ntouch[Node] = 1;
      nstack = 0;
      vtkReebGraphStackPush(Node);

      while (vtkReebGraphStackSize())
      {

        int N = vtkReebGraphStackTop();
        vtkReebGraphStackPop();

        for (int dir = 0; dir <= 1; dir++)
        {
          for (vtkIdType A = (!dir) ? (this->GetNode(N)->ArcDownId) : (this->GetNode(N)->ArcUpId);
               A; A = (!dir) ? (this->GetArc(A)->ArcDwId1) : (this->GetArc(A)->ArcDwId0))
          {
            int M = (!dir) ? (this->GetArc(A)->NodeId0) : (this->GetArc(A)->NodeId1);

            if (Atouch[A])
              continue;

            if (!Ntouch[M])
            {
              vtkReebGraphStackPush(M);
            }
            else
            {
              this->LoopNumber++;
              this->ArcLoopTable =
                (vtkIdType*)realloc(this->ArcLoopTable, sizeof(vtkIdType) * this->LoopNumber);
              this->ArcLoopTable[this->LoopNumber - 1] = A;
            }

            Atouch[A] = 1;
            Ntouch[M] = 1;
          }
        }
      }
    }
  }

  free(stack);
  free(Ntouch);
  free(Atouch);
}

//------------------------------------------------------------------------------
vtkIdType vtkReebGraph::Implementation::AddMeshVertex(vtkIdType vertexId, double scalar)
{
  static bool firstVertex = true;

  ScalarField[vertexId] = scalar;

  vtkIdType N0;
  ResizeMainNodeTable(1);

  // create a new node in the graph...
  N0 = this->MainNodeTable.FreeZone;
  this->MainNodeTable.FreeZone = this->GetNode(N0)->ArcDownId;
  ++(this->MainNodeTable.Number);
  memset(this->GetNode(N0), 0, sizeof(vtkReebNode));

  vtkReebNode* node = this->GetNode(N0);
  node->VertexId = vertexId;
  node->Value = scalar;
  node->ArcDownId = 0;
  node->ArcUpId = 0;
  node->IsFinalized = false;

  if (firstVertex)
  {
    this->MinimumScalarValue = node->Value;
    this->MaximumScalarValue = node->Value;
  }
  else
  {
    if (node->Value > this->MaximumScalarValue)
      this->MaximumScalarValue = node->Value;
    if (node->Value < this->MinimumScalarValue)
      this->MinimumScalarValue = node->Value;
  }
  firstVertex = false;

  return N0;
}

//------------------------------------------------------------------------------
vtkIdType vtkReebGraph::Implementation::FindDwLabel(vtkIdType nodeId, vtkReebLabelTag label)
{
  for (vtkIdType arcId = this->GetNode(nodeId)->ArcDownId; arcId;
       arcId = this->GetArc(arcId)->ArcDwId1)
  {
    for (vtkIdType labelId = this->GetArc(arcId)->LabelId0; labelId;
         labelId = this->GetLabel(labelId)->HNext)
    {
      if (this->GetLabel(labelId)->label == label)
        return labelId;
    }
  }
  return 0;
}

//------------------------------------------------------------------------------
vtkIdType vtkReebGraph::Implementation::FindUpLabel(vtkIdType nodeId, vtkReebLabelTag label)
{
  for (vtkIdType arcId = this->GetNode(nodeId)->ArcUpId; arcId;
       arcId = this->GetArc(arcId)->ArcDwId0)
  {
    for (vtkIdType labelId = this->GetArc(arcId)->LabelId0; labelId;
         labelId = this->GetLabel(labelId)->HNext)
    {
      if (this->GetLabel(labelId)->label == label)
        return labelId;
    }
  }
  return 0;
}

//------------------------------------------------------------------------------
void vtkReebGraph::Implementation::ResizeMainArcTable(int newSize)
{
  int oldsize, i;
  if ((this->MainArcTable.Size - this->MainArcTable.Number) < newSize)
  {
    oldsize = this->MainArcTable.Size;
    if (!this->MainArcTable.Size)
      this->MainArcTable.Size = newSize;
    while ((this->MainArcTable.Size - this->MainArcTable.Number) < newSize)
      this->MainArcTable.Size <<= 1;

    this->MainArcTable.Buffer =
      (vtkReebArc*)realloc(this->MainArcTable.Buffer, sizeof(vtkReebArc) * this->MainArcTable.Size);
    for (i = oldsize; i < this->MainArcTable.Size - 1; i++)
    {
      this->GetArc(i)->LabelId0 = i + 1;
      // clear arc
      this->GetArc(i)->LabelId1 = ((int)-2);
    }

    this->GetArc(i)->LabelId0 = this->MainArcTable.FreeZone;
    // clear arc
    this->GetArc(i)->LabelId1 = ((int)-2);
    this->MainArcTable.FreeZone = oldsize;
  }
}

//------------------------------------------------------------------------------
void vtkReebGraph::Implementation::ResizeMainLabelTable(int newSize)
{
  int oldsize, i;
  if ((this->MainLabelTable.Size - this->MainLabelTable.Number) < newSize)
  {
    oldsize = this->MainLabelTable.Size;
    if (!this->MainLabelTable.Size)
      this->MainLabelTable.Size = newSize;
    while ((this->MainLabelTable.Size - this->MainLabelTable.Number) < newSize)
      this->MainLabelTable.Size <<= 1;

    this->MainLabelTable.Buffer = (vtkReebLabel*)realloc(
      this->MainLabelTable.Buffer, sizeof(vtkReebLabel) * this->MainLabelTable.Size);

    for (i = oldsize; i < this->MainLabelTable.Size - 1; i++)
    {
      this->GetLabel(i)->ArcId = i + 1;
      this->GetLabel(i)->HNext = ((int)-2);
    }

    this->GetLabel(i)->ArcId = this->MainLabelTable.FreeZone;
    this->GetLabel(i)->HNext = ((int)-2);
    this->MainLabelTable.FreeZone = oldsize;
  }
}

//------------------------------------------------------------------------------
vtkIdType vtkReebGraph::Implementation::AddPath(
  int nodeNumber, vtkIdType* nodeOffset, vtkReebLabelTag label)
{
  vtkIdType i, Lprev, Ret = 0;

  this->ResizeMainArcTable(nodeNumber - 1);

  if (label)
    ResizeMainLabelTable(nodeNumber - 1);

  Lprev = 0;
  for (i = 0; i < (nodeNumber - 1); i++)
  {
    vtkIdType N0 = nodeOffset[i];
    vtkIdType N1 = nodeOffset[i + 1];

    // create a new arc in the graph
    int A = this->MainArcTable.FreeZone;
    this->MainArcTable.FreeZone = this->GetArc(A)->LabelId0;
    ++(this->MainArcTable.Number);
    memset(this->GetArc(A), 0, sizeof(vtkReebArc));
    vtkReebArc* arc = this->GetArc(A);

    int L = 0;

    if (!Ret)
      Ret = A;

    if (label)
    {
      vtkReebLabel* temporaryLabel;

      // create a new label in the graph
      L = this->MainLabelTable.FreeZone;
      this->MainLabelTable.FreeZone = this->GetLabel(L)->ArcId;
      ++(this->MainLabelTable.Number);
      memset(this->GetLabel(L), 0, sizeof(vtkReebLabel));
      temporaryLabel = this->GetLabel(L);

      temporaryLabel->ArcId = A;
      temporaryLabel->label = label;
      temporaryLabel->VPrev = Lprev;
    }

    arc->NodeId0 = N0;
    arc->NodeId1 = N1;
    arc->LabelId0 = arc->LabelId1 = L;

    vtkReebGraphAddUpArc(this, N0, A);
    vtkReebGraphAddDownArc(this, N1, A);

    if (label)
    {
      if (Lprev)
        this->GetLabel(Lprev)->VNext = L;
      Lprev = L;
    }
  }

  return Ret;
}

//------------------------------------------------------------------------------
void vtkReebGraph::Implementation::Collapse(vtkIdType startingNode, vtkIdType endingNode,
  vtkReebLabelTag startingLabel, vtkReebLabelTag endingLabel)
{

  int L0, L0n, L1, L1n;
  int cont[3] = { 0, 0, 0 }, Case;

  if (startingNode == endingNode)
    return;

  vtkReebNode* nstart = this->GetNode(startingNode);
  vtkReebNode* nend = this->GetNode(endingNode);

  if (!vtkReebGraphIsSmaller(this, startingNode, endingNode, nstart, nend))
  {
    vtkReebGraphSwapVars(int, startingNode, endingNode);
    vtkReebGraphSwapVars(vtkReebNode*, nstart, nend);
  }

  L0 = FindUpLabel(startingNode, startingLabel);
  L1 = FindUpLabel(startingNode, endingLabel);

  while (true)
  {
    int A0 = this->GetLabel(L0)->ArcId;
    vtkReebArc* a0 = this->GetArc(A0);
    int A1 = this->GetLabel(L1)->ArcId;
    vtkReebArc* a1 = this->GetArc(A1);

    vtkReebNode* down0 = this->GetNode(a0->NodeId0);
    vtkReebNode* up0 = this->GetNode(a0->NodeId1);
    vtkReebNode* up1 = this->GetNode(a1->NodeId1);

    /* it is the same arc, no semplification is done */
    if (A0 == A1)
    {
      Case = 0;
      L0n = this->GetLabel(L0)->VNext;
      L1n = this->GetLabel(L1)->VNext;
    }
    /* there are two arcs connecting the same start-end node */
    else if (A0 != A1 && a0->NodeId1 == a1->NodeId1)
    {
      Case = 1;
      vtkReebGraphRemoveUpArc(this, a0->NodeId0, A1);
      vtkReebGraphRemoveDownArc(this, a0->NodeId1, A1);

      for (int Lcur = this->GetArc(A1)->LabelId0; Lcur; Lcur = this->GetLabel(Lcur)->HNext)
      {
        this->GetLabel(Lcur)->ArcId = A0;
      }

      this->GetLabel(this->GetArc(A1)->LabelId0)->HPrev = this->GetArc(A0)->LabelId1;
      this->GetLabel(this->GetArc(A0)->LabelId1)->HNext = this->GetArc(A1)->LabelId0;
      this->GetArc(A0)->LabelId1 = this->GetArc(A1)->LabelId1;

      this->GetArc(A1)->LabelId0 = 0;
      this->GetArc(A1)->LabelId1 = 0;

      // delete the arc from the graph...
      this->GetArc(A1)->LabelId1 = ((int)-2);
      this->GetArc(A1)->LabelId0 = this->MainArcTable.FreeZone;
      this->MainArcTable.FreeZone = (A1);
      --(this->MainArcTable.Number);

      L0n = this->GetLabel(L0)->VNext;
      L1n = this->GetLabel(L1)->VNext;
    }
    else
    {
      if (historyOn)
      {
        vtkReebCancellation c;
        int downVertex, middleVertex, upVertex;
        downVertex = down0->VertexId;
        middleVertex = up0->VertexId;
        upVertex = up1->VertexId;
        c.removedArcs.emplace_back(downVertex, upVertex);
        c.insertedArcs.emplace_back(downVertex, middleVertex);
        c.insertedArcs.emplace_back(middleVertex, upVertex);
        this->cancellationHistory.push_back(c);
      }
      // a more complicate situation, collapse reaching the less ending point of
      // the arcs.
      Case = 2;
      {
        vtkReebNode* a0n1 = this->GetNode(a0->NodeId1);
        vtkReebNode* a1n1 = this->GetNode(a1->NodeId1);
        if (!vtkReebGraphIsSmaller(this, a0->NodeId1, a1->NodeId1, a0n1, a1n1))
        {
          vtkReebGraphSwapVars(int, A0, A1);
          vtkReebGraphSwapVars(int, L0, L1);
          vtkReebGraphSwapVars(vtkReebArc*, a0, a1);
        }
      }

      vtkReebGraphRemoveUpArc(this, a0->NodeId0, A1);
      a1->NodeId0 = a0->NodeId1;
      vtkReebGraphAddUpArc(this, a0->NodeId1, A1);

      //"replicate" labels from A1 to A0
      for (int Lcur = this->GetArc(A1)->LabelId0; Lcur; Lcur = this->GetLabel(Lcur)->HNext)
      {
        int Lnew;
        ResizeMainLabelTable(1);

        // create a new label in the graph
        Lnew = this->MainLabelTable.FreeZone;
        this->MainLabelTable.FreeZone = this->GetLabel(Lnew)->ArcId;
        ++(this->MainLabelTable.Number);
        memset(this->GetLabel(Lnew), 0, sizeof(vtkReebLabel));

        vtkReebLabel* lnew = this->GetLabel(Lnew);
        vtkReebLabel* lcur = this->GetLabel(Lcur);
        lnew->ArcId = A0;
        lnew->VPrev = lcur->VPrev;

        if (lcur->VPrev)
          this->GetLabel(lcur->VPrev)->VNext = Lnew;

        lcur->VPrev = Lnew;
        lnew->VNext = Lcur;
        lnew->label = lcur->label;

        lnew->HNext = 0;
        lnew->HPrev = this->GetArc(A0)->LabelId1;
        this->GetLabel(this->GetArc(A0)->LabelId1)->HNext = Lnew;

        this->GetArc(A0)->LabelId1 = Lnew;
      }

      L0n = this->GetLabel(L0)->VNext;
      L1n = L1;
    }

    ++cont[Case];

    int N0 = a0->NodeId0;
    vtkReebNode* n0 = this->GetNode(N0);

    if (n0->IsFinalized && vtkReebGraphIsRegular(this, n0))
    {
      if (historyOn)
      {
        vtkReebArc *up = this->GetArc(n0->ArcUpId), *down = this->GetArc(n0->ArcDownId);

        vtkReebCancellation c;
        int v0, v1, v2, v3;

        v0 = this->GetNode(up->NodeId0)->VertexId;
        v1 = this->GetNode(up->NodeId1)->VertexId;
        v2 = this->GetNode(down->NodeId0)->VertexId;
        v3 = this->GetNode(down->NodeId1)->VertexId;

        c.removedArcs.emplace_back(v0, v1);
        c.removedArcs.emplace_back(v2, v3);
        c.insertedArcs.emplace_back(v2, v1);
        this->cancellationHistory.push_back(c);
      }
      this->CollapseVertex(N0, n0);
    }

    /* end condition */
    if (a0->NodeId1 == endingNode)
    {
      vtkReebNode* nendNode = this->GetNode(endingNode);

      if (nendNode->IsFinalized && vtkReebGraphIsRegular(this, nendNode))
      {
        if (historyOn)
        {
          vtkReebArc *up = this->GetArc(this->GetNode(endingNode)->ArcUpId),
                     *down = this->GetArc(this->GetNode(endingNode)->ArcDownId);

          vtkReebCancellation c;

          int v0, v1, v2, v3;
          v0 = this->GetNode(up->NodeId0)->VertexId;
          v1 = this->GetNode(up->NodeId1)->VertexId;
          v2 = this->GetNode(down->NodeId0)->VertexId;
          v3 = this->GetNode(down->NodeId1)->VertexId;

          c.removedArcs.emplace_back(v0, v1);
          c.removedArcs.emplace_back(v2, v3);
          c.insertedArcs.emplace_back(v2, v1);
          this->cancellationHistory.push_back(c);
        }
        this->CollapseVertex(endingNode, nendNode);
      }

      return;
    }

    L0 = L0n;
    L1 = L1n;
  }
}

void vtkReebGraph::Implementation::SimplifyLabels(
  const vtkIdType nodeId, vtkReebLabelTag onlyLabel, bool goDown, bool goUp)
{
  int L, Lnext;
  vtkReebLabel* l;
  vtkReebNode* n = this->GetNode(nodeId);

  // I remove all Labels (paths) which start from me
  if (goDown)
  {
    int Anext;
    for (vtkIdType A = n->ArcDownId; A; A = Anext)
    {
      Anext = this->GetArc(A)->ArcDwId1;
      for (L = this->GetArc(A)->LabelId0; L; L = Lnext)
      {
        Lnext = this->GetLabel(L)->HNext;

        if (!(l = this->GetLabel(L))->VNext) //...starts from me!
        {
          if (!onlyLabel || onlyLabel == this->GetLabel(L)->label)
          {
            int Lprev;
            for (int Lcur = L; Lcur; Lcur = Lprev)
            {
              vtkReebLabel* lcur = this->GetLabel(Lcur);
              Lprev = lcur->VPrev;
              int CurA = lcur->ArcId;
              if (lcur->HPrev)
                this->GetLabel(lcur->HPrev)->HNext = lcur->HNext;
              else
                this->GetArc(CurA)->LabelId0 = lcur->HNext;
              if (lcur->HNext)
                this->GetLabel(lcur->HNext)->HPrev = lcur->HPrev;
              else
                this->GetArc(CurA)->LabelId1 = lcur->HPrev;

              // delete the label
              this->GetLabel(Lcur)->HNext = ((int)-2);
              this->GetLabel(Lcur)->ArcId = this->MainLabelTable.FreeZone;
              this->MainLabelTable.FreeZone = (Lcur);
              --(this->MainLabelTable.Number);
            }
          }
        }
      }
    }
  }

  // Remove all Labels (paths) which start from here

  if (goUp && !(this->GetNode(nodeId)->ArcUpId == ((int)-2)))
  {
    int Anext;
    for (vtkIdType A = n->ArcUpId; A; A = Anext)
    {
      Anext = this->GetArc(A)->ArcDwId0;
      for (L = this->GetArc(A)->LabelId0; L; L = Lnext)
      {
        Lnext = this->GetLabel(L)->HNext;

        if (!(l = this->GetLabel(L))->VPrev) //...starts from me!
        {
          if (!onlyLabel || onlyLabel == this->GetLabel(L)->label)
          {
            int myLnext;
            for (int Lcur = L; Lcur; Lcur = myLnext)
            {
              vtkReebLabel* lcur = this->GetLabel(Lcur);
              myLnext = lcur->VNext;
              int CurA = lcur->ArcId;
              vtkReebArc* cura = this->GetArc(CurA);
              if (lcur->HPrev)
              {
                this->GetLabel(lcur->HPrev)->HNext = lcur->HNext;
              }
              else
              {
                cura->LabelId0 = lcur->HNext;
              }
              if (lcur->HNext)
              {
                this->GetLabel(lcur->HNext)->HPrev = lcur->HPrev;
              }
              else
              {
                cura->LabelId1 = lcur->HPrev;
              }

              // delete the label...
              this->GetLabel(Lcur)->HNext = ((int)-2);
              this->GetLabel(Lcur)->ArcId = this->MainLabelTable.FreeZone;
              this->MainLabelTable.FreeZone = (Lcur);
              --(this->MainLabelTable.Number);
            }
          }
        }
      }
    }
  }
}

//------------------------------------------------------------------------------
void vtkReebGraph::Implementation::EndVertex(const vtkIdType N)
{
  vtkReebNode* n = this->GetNode(N);

  n->IsFinalized = true;

  if (!(this->GetNode(N)->ArcUpId == ((int)-2)))
  {
    this->SimplifyLabels(N);

    if (!(this->GetNode(N)->ArcUpId == ((int)-2)))
    {
      // special case for regular point. A node is regular if it has one
      // arc down and one arc up. In this case it can disappear

      if (vtkReebGraphIsRegular(this, n))
      {
        this->CollapseVertex(N, n);
      }
    }
  }
}

//------------------------------------------------------------------------------
int vtkReebGraph::Implementation::AddMeshTetrahedron(vtkIdType vertex0Id, double f0,
  vtkIdType vertex1Id, double f1, vtkIdType vertex2Id, double f2, vtkIdType vertex3Id, double f3)
{
  vtkIdType vertex0, vertex1, vertex2, vertex3;

  vertex0 = this->VertexStream[vertex0Id];
  vertex1 = this->VertexStream[vertex1Id];
  vertex2 = this->VertexStream[vertex2Id];
  vertex3 = this->VertexStream[vertex3Id];

  int N0 = this->VertexMap[vertex0];
  int N1 = this->VertexMap[vertex1];
  int N2 = this->VertexMap[vertex2];
  int N3 = this->VertexMap[vertex3];

  // Consistency less check
  if (f3 < f2 || (f3 == f2 && vertex3 < vertex2))
  {
    vtkReebGraphSwapVars(int, vertex2, vertex3);
    vtkReebGraphSwapVars(int, N2, N3);
    vtkReebGraphSwapVars(double, f2, f3);
  }
  if (f2 < f1 || (f2 == f1 && vertex2 < vertex1))
  {
    vtkReebGraphSwapVars(int, vertex1, vertex2);
    vtkReebGraphSwapVars(int, N1, N2);
    vtkReebGraphSwapVars(double, f1, f2);
  }
  if (f1 < f0 || (f1 == f0 && vertex1 < vertex0))
  {
    vtkReebGraphSwapVars(int, vertex0, vertex1);
    vtkReebGraphSwapVars(int, N0, N1);
    vtkReebGraphSwapVars(double, f0, f1);
  }
  if (f3 < f2 || (f3 == f2 && vertex3 < vertex2))
  {
    vtkReebGraphSwapVars(int, vertex2, vertex3);
    vtkReebGraphSwapVars(int, N2, N3);
    vtkReebGraphSwapVars(double, f2, f3);
  }
  if (f2 < f1 || (f2 == f1 && vertex2 < vertex1))
  {
    vtkReebGraphSwapVars(int, vertex1, vertex2);
    vtkReebGraphSwapVars(int, N1, N2);
    vtkReebGraphSwapVars(double, f1, f2);
  }
  if (f3 < f2 || (f3 == f2 && vertex3 < vertex2))
  {
    vtkReebGraphSwapVars(int, vertex2, vertex3);
    vtkReebGraphSwapVars(int, N2, N3);
    vtkReebGraphSwapVars(double, f2, f3);
  }

  vtkIdType t0[] = { vertex0, vertex1, vertex2 }, t1[] = { vertex0, vertex1, vertex3 },
            t2[] = { vertex0, vertex2, vertex3 }, t3[] = { vertex1, vertex2, vertex3 };
  vtkIdType* cellIds[4];
  cellIds[0] = t0;
  cellIds[1] = t1;
  cellIds[2] = t2;
  cellIds[3] = t3;

  for (int i = 0; i < 3; i++)
  {
    int n0 = this->VertexMap[cellIds[i][0]], n1 = this->VertexMap[cellIds[i][1]],
        n2 = this->VertexMap[cellIds[i][2]];

    vtkReebLabelTag Label01 =
      ((vtkReebLabelTag)cellIds[i][0]) | (((vtkReebLabelTag)cellIds[i][1]) << 32);
    vtkReebLabelTag Label12 =
      ((vtkReebLabelTag)cellIds[i][1]) | (((vtkReebLabelTag)cellIds[i][2]) << 32);
    vtkReebLabelTag Label02 =
      ((vtkReebLabelTag)cellIds[i][0]) | (((vtkReebLabelTag)cellIds[i][2]) << 32);

    if (!this->FindUpLabel(n0, Label01))
    {
      vtkIdType N01[] = { n0, n1 };
      this->AddPath(2, N01, Label01);
    }
    if (!this->FindUpLabel(n1, Label12))
    {
      vtkIdType N12[] = { n1, n2 };
      this->AddPath(2, N12, Label12);
    }
    if (!this->FindUpLabel(n0, Label02))
    {
      vtkIdType N02[] = { n0, n2 };
      this->AddPath(2, N02, Label02);
    }

    this->Collapse(n0, n1, Label01, Label02);
    this->Collapse(n1, n2, Label12, Label02);
  }

  if (!(--(this->TriangleVertexMap[vertex0])))
    this->EndVertex(N0);
  if (!(--(this->TriangleVertexMap[vertex1])))
    this->EndVertex(N1);
  if (!(--(this->TriangleVertexMap[vertex2])))
    this->EndVertex(N2);
  if (!(--(this->TriangleVertexMap[vertex3])))
    this->EndVertex(N3);

  return 1;
}

//------------------------------------------------------------------------------
int vtkReebGraph::Implementation::AddMeshTriangle(
  vtkIdType vertex0Id, double f0, vtkIdType vertex1Id, double f1, vtkIdType vertex2Id, double f2)
{

  int vertex0 = this->VertexStream[vertex0Id], vertex1 = this->VertexStream[vertex1Id],
      vertex2 = this->VertexStream[vertex2Id];

  int N0 = this->VertexMap[vertex0];
  int N1 = this->VertexMap[vertex1];
  int N2 = this->VertexMap[vertex2];

  // Consistency less check
  if (f2 < f1 || (f2 == f1 && vertex2 < vertex1))
  {
    vtkReebGraphSwapVars(int, vertex1, vertex2);
    vtkReebGraphSwapVars(int, N1, N2);
    vtkReebGraphSwapVars(double, f1, f2);
  }
  if (f1 < f0 || (f1 == f0 && vertex1 < vertex0))
  {
    vtkReebGraphSwapVars(int, vertex0, vertex1);
    vtkReebGraphSwapVars(int, N0, N1);
    // there is a useless assignment to f0 below
    // keeping for clarity and consistency
    vtkReebGraphSwapVars(double, f0, f1);
  }
  if (f2 < f1 || (f2 == f1 && vertex2 < vertex1))
  {
    vtkReebGraphSwapVars(int, vertex1, vertex2);
    vtkReebGraphSwapVars(int, N1, N2);
    // there is a useless assignment below
    // keeping for clarity and consistency
    vtkReebGraphSwapVars(double, f1, f2);
  }

  vtkReebLabelTag Label01 = ((vtkReebLabelTag)vertex0) | (((vtkReebLabelTag)vertex1) << 32);
  vtkReebLabelTag Label12 = ((vtkReebLabelTag)vertex1) | (((vtkReebLabelTag)vertex2) << 32);
  vtkReebLabelTag Label02 = ((vtkReebLabelTag)vertex0) | (((vtkReebLabelTag)vertex2) << 32);

  if (!this->FindUpLabel(N0, Label01))
  {
    vtkIdType N01[] = { N0, N1 };
    this->AddPath(2, N01, Label01);
  }
  if (!this->FindUpLabel(N1, Label12))
  {
    vtkIdType N12[] = { N1, N2 };
    this->AddPath(2, N12, Label12);
  }
  if (!this->FindUpLabel(N0, Label02))
  {
    vtkIdType N02[] = { N0, N2 };
    this->AddPath(2, N02, Label02);
  }

  this->Collapse(N0, N1, Label01, Label02);
  this->Collapse(N1, N2, Label12, Label02);

  if (!(--(this->TriangleVertexMap[vertex0])))
    this->EndVertex(N0);
  if (!(--(this->TriangleVertexMap[vertex1])))
    this->EndVertex(N1);
  if (!(--(this->TriangleVertexMap[vertex2])))
    this->EndVertex(N2);

  return 1;
}

//------------------------------------------------------------------------------
int vtkReebGraph::StreamTetrahedron(vtkIdType vertex0Id, double scalar0, vtkIdType vertex1Id,
  double scalar1, vtkIdType vertex2Id, double scalar2, vtkIdType vertex3Id, double scalar3)
{
  return this->Storage->StreamTetrahedron(
    vertex0Id, scalar0, vertex1Id, scalar1, vertex2Id, scalar2, vertex3Id, scalar3);
}

//------------------------------------------------------------------------------
int vtkReebGraph::Implementation::StreamTetrahedron(vtkIdType vertex0Id, double scalar0,
  vtkIdType vertex1Id, double scalar1, vtkIdType vertex2Id, double scalar2, vtkIdType vertex3Id,
  double scalar3)
{
  if (!this->VertexMapAllocatedSize)
  {
    // first allocate an arbitrary size
    this->VertexMapAllocatedSize = vtkReebGraphInitialStreamSize;
    this->VertexMap = (vtkIdType*)malloc(sizeof(vtkIdType) * this->VertexMapAllocatedSize);
    memset(this->VertexMap, 0, sizeof(vtkIdType) * this->VertexMapAllocatedSize);
    this->VertexStream.clear();
  }
  else if (this->VertexMapSize >= this->VertexMapAllocatedSize - 4)
  {
    int oldSize = this->VertexMapAllocatedSize;
    this->VertexMapAllocatedSize <<= 1;
    this->VertexMap =
      (vtkIdType*)realloc(this->VertexMap, sizeof(vtkIdType) * this->VertexMapAllocatedSize);
    for (int i = oldSize; i < this->VertexMapAllocatedSize - 1; i++)
      this->VertexMap[i] = 0;
  }

  // same thing with the triangle map
  if (!this->TriangleVertexMapAllocatedSize)
  {
    // first allocate an arbitrary size
    this->TriangleVertexMapAllocatedSize = vtkReebGraphInitialStreamSize;
    this->TriangleVertexMap = (int*)malloc(sizeof(int) * this->TriangleVertexMapAllocatedSize);
    memset(this->TriangleVertexMap, 0, sizeof(int) * this->TriangleVertexMapAllocatedSize);
  }
  else if (this->TriangleVertexMapSize >= this->TriangleVertexMapAllocatedSize - 4)
  {
    int oldSize = this->TriangleVertexMapAllocatedSize;
    this->TriangleVertexMapAllocatedSize <<= 1;
    this->TriangleVertexMap =
      (int*)realloc(this->TriangleVertexMap, sizeof(int) * this->TriangleVertexMapAllocatedSize);
    for (int i = oldSize; i < this->TriangleVertexMapAllocatedSize - 1; i++)
      this->TriangleVertexMap[i] = 0;
  }

  // Add the vertices to the stream
  std::map<int, int>::iterator sIter;

  // vertex0
  sIter = this->VertexStream.find(vertex0Id);
  if (sIter == this->VertexStream.end())
  {
    // this vertex hasn't been streamed yet, let's add it
    this->VertexStream[vertex0Id] = this->VertexMapSize;
    this->VertexMap[this->VertexMapSize] = this->AddMeshVertex(vertex0Id, scalar0);
    this->VertexMapSize++;
    this->TriangleVertexMapSize++;
  }

  // vertex1
  sIter = this->VertexStream.find(vertex1Id);
  if (sIter == this->VertexStream.end())
  {
    // this vertex hasn't been streamed yet, let's add it
    this->VertexStream[vertex1Id] = this->VertexMapSize;
    this->VertexMap[this->VertexMapSize] = this->AddMeshVertex(vertex1Id, scalar1);
    this->VertexMapSize++;
    this->TriangleVertexMapSize++;
  }

  // vertex2
  sIter = this->VertexStream.find(vertex2Id);
  if (sIter == this->VertexStream.end())
  {
    // this vertex hasn't been streamed yet, let's add it
    this->VertexStream[vertex2Id] = this->VertexMapSize;
    this->VertexMap[this->VertexMapSize] = this->AddMeshVertex(vertex2Id, scalar2);
    this->VertexMapSize++;
    this->TriangleVertexMapSize++;
  }

  // vertex3
  sIter = this->VertexStream.find(vertex3Id);
  if (sIter == this->VertexStream.end())
  {
    // this vertex hasn't been streamed yet, let's add it
    this->VertexStream[vertex3Id] = this->VertexMapSize;
    this->VertexMap[this->VertexMapSize] = this->AddMeshVertex(vertex3Id, scalar3);
    this->VertexMapSize++;
    this->TriangleVertexMapSize++;
  }

  this->AddMeshTetrahedron(
    vertex0Id, scalar0, vertex1Id, scalar1, vertex2Id, scalar2, vertex3Id, scalar3);

  return 0;
}

//------------------------------------------------------------------------------
int vtkReebGraph::StreamTriangle(vtkIdType vertex0Id, double scalar0, vtkIdType vertex1Id,
  double scalar1, vtkIdType vertex2Id, double scalar2)
{
  return this->Storage->StreamTriangle(vertex0Id, scalar0, vertex1Id, scalar1, vertex2Id, scalar2);
}

//------------------------------------------------------------------------------
int vtkReebGraph::Implementation::StreamTriangle(vtkIdType vertex0Id, double scalar0,
  vtkIdType vertex1Id, double scalar1, vtkIdType vertex2Id, double scalar2)
{
  if (!this->VertexMapAllocatedSize)
  {
    // first allocate an arbitrary size
    this->VertexMapAllocatedSize = vtkReebGraphInitialStreamSize;
    this->VertexMap = (vtkIdType*)malloc(sizeof(vtkIdType) * this->VertexMapAllocatedSize);
    memset(this->VertexMap, 0, sizeof(vtkIdType) * this->VertexMapAllocatedSize);
  }
  else if (this->VertexMapSize >= this->VertexMapAllocatedSize - 3)
  {
    int oldSize = this->VertexMapAllocatedSize;
    this->VertexMapAllocatedSize <<= 1;
    this->VertexMap =
      (vtkIdType*)realloc(this->VertexMap, sizeof(vtkIdType) * this->VertexMapAllocatedSize);
    for (int i = oldSize; i < this->VertexMapAllocatedSize - 1; i++)
      this->VertexMap[i] = 0;
  }

  // same thing with the triangle map
  if (!this->TriangleVertexMapAllocatedSize)
  {
    // first allocate an arbitrary size
    this->TriangleVertexMapAllocatedSize = vtkReebGraphInitialStreamSize;
    this->TriangleVertexMap = (int*)malloc(sizeof(int) * this->TriangleVertexMapAllocatedSize);
    memset(this->TriangleVertexMap, 0, sizeof(int) * this->TriangleVertexMapAllocatedSize);
  }
  else if (this->TriangleVertexMapSize >= this->TriangleVertexMapAllocatedSize - 3)
  {
    int oldSize = this->TriangleVertexMapAllocatedSize;
    this->TriangleVertexMapAllocatedSize <<= 1;
    this->TriangleVertexMap =
      (int*)realloc(this->TriangleVertexMap, sizeof(int) * this->TriangleVertexMapAllocatedSize);
    for (int i = oldSize; i < this->TriangleVertexMapAllocatedSize - 1; i++)
      this->TriangleVertexMap[i] = 0;
  }

  // Add the vertices to the stream
  std::map<int, int>::iterator sIter;

  // vertex0
  sIter = this->VertexStream.find(vertex0Id);
  if (sIter == this->VertexStream.end())
  {
    // this vertex hasn't been streamed yet, let's add it
    this->VertexStream[vertex0Id] = this->VertexMapSize;
    this->VertexMap[this->VertexMapSize] = this->AddMeshVertex(vertex0Id, scalar0);
    this->VertexMapSize++;
    this->TriangleVertexMapSize++;
  }

  // vertex1
  sIter = this->VertexStream.find(vertex1Id);
  if (sIter == this->VertexStream.end())
  {
    // this vertex hasn't been streamed yet, let's add it
    this->VertexStream[vertex1Id] = this->VertexMapSize;
    this->VertexMap[this->VertexMapSize] = this->AddMeshVertex(vertex1Id, scalar1);
    this->VertexMapSize++;
    this->TriangleVertexMapSize++;
  }

  // vertex2
  sIter = this->VertexStream.find(vertex2Id);
  if (sIter == this->VertexStream.end())
  {
    // this vertex hasn't been streamed yet, let's add it
    this->VertexStream[vertex2Id] = this->VertexMapSize;
    this->VertexMap[this->VertexMapSize] = this->AddMeshVertex(vertex2Id, scalar2);
    this->VertexMapSize++;
    this->TriangleVertexMapSize++;
  }

  this->AddMeshTriangle(vertex0Id, scalar0, vertex1Id, scalar1, vertex2Id, scalar2);

  return 0;
}

//------------------------------------------------------------------------------
int vtkReebGraph::Build(vtkPolyData* mesh, vtkDataArray* scalarField)
{
  for (vtkIdType i = 0; i < mesh->GetNumberOfCells(); i++)
  {
    vtkCell* triangle = mesh->GetCell(i);
    vtkIdList* trianglePointList = triangle->GetPointIds();
    if (trianglePointList->GetNumberOfIds() != 3)
      return vtkReebGraph::ERR_NOT_A_SIMPLICIAL_MESH;
    StreamTriangle(trianglePointList->GetId(0),
      scalarField->GetComponent(trianglePointList->GetId(0), 0), trianglePointList->GetId(1),
      scalarField->GetComponent(trianglePointList->GetId(1), 0), trianglePointList->GetId(2),
      scalarField->GetComponent(trianglePointList->GetId(2), 0));
  }

  this->Storage->inputMesh = mesh;
  this->Storage->inputScalarField = scalarField;

  this->CloseStream();

  return 0;
}

//------------------------------------------------------------------------------
int vtkReebGraph::Build(vtkUnstructuredGrid* mesh, vtkDataArray* scalarField)
{
  for (vtkIdType i = 0; i < mesh->GetNumberOfCells(); i++)
  {
    vtkCell* tet = mesh->GetCell(i);
    vtkIdList* tetPointList = tet->GetPointIds();
    if (tetPointList->GetNumberOfIds() != 4)
      return vtkReebGraph::ERR_NOT_A_SIMPLICIAL_MESH;
    StreamTetrahedron(tetPointList->GetId(0), scalarField->GetComponent(tetPointList->GetId(0), 0),
      tetPointList->GetId(1), scalarField->GetComponent(tetPointList->GetId(1), 0),
      tetPointList->GetId(2), scalarField->GetComponent(tetPointList->GetId(2), 0),
      tetPointList->GetId(3), scalarField->GetComponent(tetPointList->GetId(3), 0));
  }

  this->Storage->inputMesh = mesh;
  this->Storage->inputScalarField = scalarField;

  this->CloseStream();

  return 0;
}

//------------------------------------------------------------------------------
int vtkReebGraph::Implementation::GetNumberOfArcs()
{
  if (!this->ArcNumber)
    for (vtkIdType arcId = 1; arcId < this->MainArcTable.Size; arcId++)
    {
      // check if arc is cleared
      if (!(this->GetArc(arcId)->LabelId1 == ((int)-2)))
        this->ArcNumber++;
    }

  return this->ArcNumber;
}

//------------------------------------------------------------------------------
int vtkReebGraph::Implementation::GetNumberOfConnectedComponents()
{
  if (!this->ArcLoopTable)
    this->FindLoops();
  return this->ConnectedComponentNumber;
}

//------------------------------------------------------------------------------
int vtkReebGraph::Implementation::GetNumberOfNodes()
{
  if (!this->NodeNumber)
    for (vtkIdType nodeId = 1; nodeId < this->MainNodeTable.Size; nodeId++)
    {
      // check if node is cleared
      if (!(this->GetNode(nodeId)->ArcUpId == ((int)-2)))
        this->NodeNumber++;
    }

  return this->NodeNumber;
}

//------------------------------------------------------------------------------
vtkIdType vtkReebGraph::Implementation::GetNextNodeId()
{
  for (vtkIdType nodeId = this->currentNodeId + 1; nodeId < this->MainNodeTable.Size; nodeId++)
  {
    // check if node is cleared
    if (!(this->GetNode(nodeId)->ArcUpId == ((int)-2)))
    {
      this->currentNodeId = nodeId;
      return this->currentNodeId;
    }
  }

  return this->currentNodeId;
}

//------------------------------------------------------------------------------
vtkIdType vtkReebGraph::Implementation::GetPreviousNodeId()
{
  if (!this->currentNodeId)
  {
    return this->GetNextNodeId();
  }

  for (vtkIdType nodeId = this->currentNodeId - 1; nodeId > 0; nodeId--)
  {
    // check if node is cleared
    if (!(this->GetNode(nodeId)->ArcUpId == ((int)-2)))
    {
      this->currentNodeId = nodeId;
      return this->currentNodeId;
    }
  }

  return this->currentNodeId;
}

//------------------------------------------------------------------------------
vtkIdType vtkReebGraph::Implementation::GetNextArcId()
{
  for (vtkIdType arcId = this->currentArcId + 1; arcId < this->MainArcTable.Size; arcId++)
  {
    // check if arc is cleared
    if (!(this->GetArc(arcId)->LabelId1 == ((int)-2)))
    {
      this->currentArcId = arcId;
      return this->currentArcId;
    }
  }

  return this->currentArcId;
}

//------------------------------------------------------------------------------
vtkIdType vtkReebGraph::Implementation::GetPreviousArcId()
{
  if (!this->currentArcId)
  {
    return this->GetNextArcId();
  }

  for (vtkIdType arcId = this->currentArcId - 1; arcId > 0; arcId--)
  {
    // check if arc is cleared
    if (!(this->GetArc(arcId)->LabelId1 == ((int)-2)))
    {
      this->currentArcId = arcId;
      return this->currentArcId;
    }
  }

  return this->currentArcId;
}

//------------------------------------------------------------------------------
vtkIdType vtkReebGraph::Implementation::GetArcDownNodeId(vtkIdType arcId)
{
  return (this->GetArc(arcId))->NodeId0;
}

//------------------------------------------------------------------------------
vtkIdType vtkReebGraph::Implementation::GetArcUpNodeId(vtkIdType arcId)
{
  return (this->GetArc(arcId))->NodeId1;
}

//------------------------------------------------------------------------------
double vtkReebGraph::Implementation::GetNodeScalarValue(vtkIdType nodeId)
{
  return (this->GetNode(nodeId))->Value;
}

//------------------------------------------------------------------------------
vtkIdType vtkReebGraph::Implementation::GetNodeVertexId(vtkIdType nodeId)
{
  return (this->GetNode(nodeId))->VertexId;
}

//------------------------------------------------------------------------------
int vtkReebGraph::Build(vtkPolyData* mesh, vtkIdType scalarFieldId)
{
  vtkPointData* pointData = mesh->GetPointData();
  vtkDataArray* scalarField = pointData->GetArray(scalarFieldId);

  if (!scalarField)
    return vtkReebGraph::ERR_NO_SUCH_FIELD;

  return this->Build(mesh, scalarField);
}

//------------------------------------------------------------------------------
int vtkReebGraph::Build(vtkUnstructuredGrid* mesh, vtkIdType scalarFieldId)
{
  vtkPointData* pointData = mesh->GetPointData();
  vtkDataArray* scalarField = pointData->GetArray(scalarFieldId);

  if (!scalarField)
    return vtkReebGraph::ERR_NO_SUCH_FIELD;

  return this->Build(mesh, scalarField);
}

//------------------------------------------------------------------------------
int vtkReebGraph::Build(vtkPolyData* mesh, const char* scalarFieldName)
{
  int scalarFieldId = 0;

  vtkPointData* pointData = mesh->GetPointData();
  vtkDataArray* scalarField = pointData->GetArray(scalarFieldName, scalarFieldId);

  if (!scalarField)
    return vtkReebGraph::ERR_NO_SUCH_FIELD;

  return this->Build(mesh, scalarField);
}

//------------------------------------------------------------------------------
int vtkReebGraph::Build(vtkUnstructuredGrid* mesh, const char* scalarFieldName)
{
  int scalarFieldId = 0;

  vtkPointData* pointData = mesh->GetPointData();
  vtkDataArray* scalarField = pointData->GetArray(scalarFieldName, scalarFieldId);

  if (!scalarField)
    return vtkReebGraph::ERR_NO_SUCH_FIELD;

  return this->Build(mesh, scalarField);
}

//------------------------------------------------------------------------------
int vtkReebGraph::Implementation::GetNumberOfLoops()
{
  if (!this->ArcLoopTable)
    this->FindLoops();
  return this->LoopNumber - this->RemovedLoopNumber;
}

//------------------------------------------------------------------------------
inline vtkIdType vtkReebGraph::Implementation::AddArc(vtkIdType nodeId0, vtkIdType nodeId1)
{
  if (!vtkReebGraphIsSmaller(
        this, nodeId0, nodeId1, this->GetNode(nodeId0), this->GetNode(nodeId1)))
    vtkReebGraphSwapVars(vtkIdType, nodeId0, nodeId1);
  vtkIdType nodevtkReebArcble[] = { nodeId0, nodeId1 };
  return AddPath(2, nodevtkReebArcble, 0);
}