#include "trianglepatterngeometry.hpp"
#include "trianglepattterntopology.hpp"
#include <algorithm>
#include <iterator>
#include <numeric>
#include <vcg/complex/algorithms/update/position.h>
#include <vcg/simplex/edge/topology.h>
#include <vcg/space/intersection2.h>
#include <wrap/io_trimesh/export.h>
/* include the support for half edges */
#include <vcg/complex/algorithms/update/halfedge_indexed.h>

size_t PatternGeometry::computeTiledValence(
    const size_t &nodeIndex,
    const std::vector<size_t> &numberOfNodesPerSlot) const {
  std::vector<PatternGeometry::EdgeType *> connectedEdges;
  vcg::edge::VEStarVE(&vert[nodeIndex], connectedEdges);
  const size_t nodeValence = connectedEdges.size();
  assert(nodeToSlotMap.count(nodeIndex) != 0);
  const size_t nodeSlotIndex = nodeToSlotMap.at(nodeIndex);
  if (nodeSlotIndex == 0) {
    return nodeValence * fanSize;
  } else if (nodeSlotIndex == 1 || nodeSlotIndex == 2) {
    size_t correspondingNodeIndex;
    if (nodeSlotIndex == 1) {
      correspondingNodeIndex = nodeIndex + 1;
    } else {
      correspondingNodeIndex = nodeIndex - 1;
    }
    std::vector<PatternGeometry::EdgeType *> connectedEdgesCorrespondingNode;
    vcg::edge::VEStarVE(&vert[correspondingNodeIndex],
                        connectedEdgesCorrespondingNode);
    size_t correspondingNodeValence = connectedEdgesCorrespondingNode.size();
    return fanSize / 2 * nodeValence + fanSize / 2 * correspondingNodeValence;
  } else if (nodeSlotIndex == 3 || nodeSlotIndex == 5) {
    size_t correspondingNodeIndex;
    size_t numberOfNodesBefore;
    size_t numberOfNodesAfter;
    if (nodeSlotIndex == 3) {
      numberOfNodesBefore =
          nodeIndex - std::accumulate(numberOfNodesPerSlot.begin(),
                                      numberOfNodesPerSlot.begin() + 3, 0);
      correspondingNodeIndex =

          std::accumulate(numberOfNodesPerSlot.begin(),
                          numberOfNodesPerSlot.begin() + 6, 0) -
          1 - numberOfNodesBefore;
    } else {
      numberOfNodesAfter =
          std::accumulate(numberOfNodesPerSlot.begin(),
                          numberOfNodesPerSlot.begin() + 6, 0) -
          1 - nodeIndex;
      correspondingNodeIndex =
          numberOfNodesAfter + std::accumulate(numberOfNodesPerSlot.begin(),
                                               numberOfNodesPerSlot.begin() + 3,
                                               0);
    }
    assert(correspondingNodeIndex < vn);
    std::vector<PatternGeometry::EdgeType *> connectedEdgesCorrespondingNode;
    vcg::edge::VEStarVE(&vert[correspondingNodeIndex],
                        connectedEdgesCorrespondingNode);
    size_t correspondingNodeValence = connectedEdgesCorrespondingNode.size();
    return nodeValence + correspondingNodeValence;
  } else if (nodeSlotIndex == 4) {
    return 2 * nodeValence;
  } else if (nodeSlotIndex == 6) {
    return nodeValence;
  } else {
    std::cerr << "Error in slot index computation" << std::endl;
  }
  assert(false);
  return 0;
}

size_t PatternGeometry::getFanSize() const { return fanSize; }

double PatternGeometry::getTriangleEdgeSize() const { return triangleEdgeSize; }

PatternGeometry::PatternGeometry() {}

std::vector<vcg::Point3d> PatternGeometry::getVertices() const {
  std::vector<VCGEdgeMesh::CoordType> verts(VN());
  for (size_t vi = 0; vi < VN(); vi++) {
    verts[vi] = vert[vi].cP();
  }
  return verts;
}

PatternGeometry PatternGeometry::createTile(PatternGeometry &pattern)
{
    const size_t fanSize = PatternGeometry().getFanSize();
    PatternGeometry fan(createFan(pattern));
    PatternGeometry tile(fan);

    if (fanSize % 2 == 1) {
        vcg::Matrix44d R;
        auto rotationAxis = vcg::Point3d(0, 0, 1);
        R.SetRotateDeg(180, rotationAxis);
        vcg::tri::UpdatePosition<PatternGeometry>::Matrix(fan, R);
    }
    vcg::Matrix44d T;
    //    const double centerAngle = 2 * M_PI / fanSize;
    //    const double triangleHeight = std::sin((M_PI - centerAngle) / 2)
    //                                  * PatternGeometry().triangleEdgeSize;
    vcg::tri::UpdateBounding<PatternGeometry>::Box(fan);
    const double triangleHeight = fan.bbox.DimY() / 2;
    T.SetTranslate(0, -2 * triangleHeight, 0);
    vcg::tri::UpdatePosition<PatternGeometry>::Matrix(fan, T);
    //    fan.registerForDrawing();
    //    polyscope::show();

    PatternGeometry fanOfFan = createFan(fan);
    vcg::tri::Append<PatternGeometry, PatternGeometry>::Mesh(tile, fanOfFan);
    vcg::tri::Clean<PatternGeometry>::MergeCloseVertex(tile, 0.0000005);
    vcg::tri::Allocator<PatternGeometry>::CompactEveryVector(tile);
    vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(tile);
    vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(tile);

    //  for (size_t vi = 0; vi < pattern.vn; vi++) {
    //    tile.vert[vi].C() = vcg::Color4b::Blue;
    //  }

    tile.setLabel("tilled_" + pattern.getLabel());
    tile.updateEigenEdgeAndVertices();
    return tile;
}

PatternGeometry PatternGeometry::createFan(PatternGeometry &pattern) {
  const size_t fanSize = PatternGeometry().getFanSize();
  PatternGeometry fan(pattern);
  PatternGeometry rotatedPattern(pattern);
  for (int rotationCounter = 1; rotationCounter < fanSize; rotationCounter++) {
    vcg::Matrix44d R;
    auto rotationAxis = vcg::Point3d(0, 0, 1);
    R.SetRotateDeg(360 / fanSize, rotationAxis);
    vcg::tri::UpdatePosition<PatternGeometry>::Matrix(rotatedPattern, R);
    vcg::tri::Append<PatternGeometry, PatternGeometry>::Mesh(fan,
                                                             rotatedPattern);
  }
  return fan;
}

PatternGeometry::PatternGeometry(PatternGeometry &other) {
  vcg::tri::Append<PatternGeometry, PatternGeometry>::MeshCopy(*this, other);
  this->vertices = other.getVertices();
  baseTriangle = other.getBaseTriangle();
  baseTriangleHeight = computeBaseTriangleHeight();
  vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
  vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(*this);
}

bool PatternGeometry::load(const std::filesystem::__cxx11::path &meshFilePath)
{
    if (!VCGEdgeMesh::load(meshFilePath)) {
        return false;
    }
    vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
    baseTriangleHeight = computeBaseTriangleHeight();
    baseTriangle = computeBaseTriangle();

    updateEigenEdgeAndVertices();
    return true;
}

void PatternGeometry::add(const std::vector<vcg::Point3d> &vertices) {
  this->vertices = vertices;
  std::for_each(vertices.begin(), vertices.end(), [&](const vcg::Point3d &p) {
    vcg::tri::Allocator<PatternGeometry>::AddVertex(*this, p);
  });
  vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
  vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(*this);
  updateEigenEdgeAndVertices();
}

void PatternGeometry::add(const std::vector<vcg::Point2i> &edges) {
  std::for_each(edges.begin(), edges.end(), [&](const vcg::Point2i &e) {
    vcg::tri::Allocator<PatternGeometry>::AddEdge(*this, e[0], e[1]);
  });
  vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
  vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(*this);
  updateEigenEdgeAndVertices();
}

void PatternGeometry::add(const std::vector<vcg::Point3d> &vertices,
                          const std::vector<vcg::Point2i> &edges)
{
    add(vertices);
    add(edges);
    updateEigenEdgeAndVertices();
}

void PatternGeometry::add(const std::vector<size_t> &numberOfNodesPerSlot,
                          const std::vector<vcg::Point2i> &edges) {
  assert(numberOfNodesPerSlot.size() == 7);
  auto vertices =
      constructVertexVector(numberOfNodesPerSlot, fanSize, triangleEdgeSize);
  add(vertices, edges);
  vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
  vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(*this);
  updateEigenEdgeAndVertices();
}

std::vector<vcg::Point3d> PatternGeometry::constructVertexVector(
    const std::vector<size_t> &numberOfNodesPerSlot, const size_t &fanSize,
    const double &triangleEdgeSize) {

  std::vector<vcg::Point3d> vertices;
  const double centerAngle = 2 * M_PI / fanSize;
  const double triangleHeight =
      std::sin((M_PI - centerAngle) / 2) * triangleEdgeSize;
  const double baseEdgeSize =
      2 * triangleEdgeSize * std::cos((M_PI - centerAngle) / 2);
  const vcg::Point3d triangleV0 = vcg::Point3d(0, 0, 0);
  const vcg::Point3d triangleV1 =
      vcg::Point3d(-baseEdgeSize / 2, -triangleHeight, 0);
  const vcg::Point3d triangleV2 =
      vcg::Point3d(baseEdgeSize / 2, -triangleHeight, 0);

  // Nodes
  if (numberOfNodesPerSlot[0] == 1) {
    //    vertices[0] = triangleV0;
    vertices.push_back(triangleV0);
  }
  if (numberOfNodesPerSlot[1] == 1) {
    //    vertices[1] = triangleV1;
    vertices.push_back(triangleV1);
  }
  if (numberOfNodesPerSlot[2] == 1) {
    //    vertices[2] = triangleV2;
    vertices.push_back(triangleV2);
  }

  // Edges
  if (numberOfNodesPerSlot[3] != 0) {
    const double offset0 = 1.0 / (numberOfNodesPerSlot[3] + 1);
    const vcg::Point3d edgeVector0(triangleV1 - triangleV0);
    for (size_t vertexIndex = 0; vertexIndex < numberOfNodesPerSlot[3];
         vertexIndex++) {

      //      vertices[std::accumulate(numberOfNodesPerSlot.begin(),
      //                               numberOfNodesPerSlot.begin() + 2, 0) +
      //               vertexIndex] =
      vertices.push_back(triangleV0 +
                         edgeVector0.operator*((vertexIndex + 1) * offset0));
    }
  }

  if (numberOfNodesPerSlot[4] == 1) {
    vertices.push_back(vcg::Point3d(0, -triangleHeight, 0));
  } else if (numberOfNodesPerSlot[4] != 0) {
    const double offset1 = 1.0 / (numberOfNodesPerSlot[4] + 1);
    const vcg::Point3d edgeVector1(triangleV2 - triangleV1);
    for (size_t vertexIndex = 0; vertexIndex < numberOfNodesPerSlot[4];
         vertexIndex++) {

      //      vertices[std::accumulate(numberOfNodesPerSlot.begin(),
      //                               numberOfNodesPerSlot.begin() + 3, 0) +
      //               vertexIndex] =
      vertices.push_back(triangleV1 +
                         edgeVector1.operator*((vertexIndex + 1) * offset1));
    }
  }

  if (numberOfNodesPerSlot[5] != 0) {
    const double offset2 = 1.0 / (numberOfNodesPerSlot[5] + 1);
    const vcg::Point3d edgeVector2(triangleV0 - triangleV2);
    for (size_t vertexIndex = 0; vertexIndex < numberOfNodesPerSlot[5];
         vertexIndex++) {
      //      vertices[std::accumulate(numberOfNodesPerSlot.begin(),
      //                               numberOfNodesPerSlot.begin() + 4, 0) +
      //               vertexIndex] =
      vertices.push_back(triangleV2 +
                         edgeVector2.operator*((vertexIndex + 1) * offset2));
    }
  }

  // Face
  if (numberOfNodesPerSlot[6] != 0) {
    const vcg::Point3d triangleCenter((triangleV0 + triangleV1 + triangleV2) /
                                      3);
    const double radius = 0.1;
    const double offsetRad = 2 * M_PI / numberOfNodesPerSlot[6];
    for (size_t vertexIndex = 0; vertexIndex < numberOfNodesPerSlot[6];
         vertexIndex++) {
      const double pX =
          triangleCenter[0] + radius * std::cos(vertexIndex * offsetRad);
      const double pY =
          triangleCenter[1] + radius * std::sin(vertexIndex * offsetRad);
      /*vertices[std::accumulate(numberOfNodesPerSlot.begin(),
                               numberOfNodesPerSlot.begin() + 5, 0) +
               vertexIndex] =*/
      vertices.push_back(vcg::Point3d(pX, pY, 0));
    }
  }
  return vertices;
}

void PatternGeometry::constructNodeToTiledValenceMap(
    const std::vector<size_t> &numberOfNodesPerSlot) {
  for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) {
    const size_t tiledValence =
        computeTiledValence(nodeIndex, numberOfNodesPerSlot);
    nodeTiledValence[nodeIndex] = tiledValence;
  }
}

std::vector<VectorType> PatternGeometry::getEdgeVectorsWithVertexAsOrigin(
    std::vector<EdgePointer> &edgePointers, const int &vi)
{
    std::vector<VectorType> incidentElementsVectors(edgePointers.size());
    for (int incidentElementsIndex = 0; incidentElementsIndex < edgePointers.size();
         incidentElementsIndex++) {
        assert(vi == getIndex(edgePointers[incidentElementsIndex]->cV(0))
               || vi == getIndex(edgePointers[incidentElementsIndex]->cV(1)));

        incidentElementsVectors[incidentElementsIndex]
            = vi == getIndex(edgePointers[incidentElementsIndex]->cV(0))
                  ? edgePointers[incidentElementsIndex]->cP(1)
                        - edgePointers[incidentElementsIndex]->cP(0)
                  : edgePointers[incidentElementsIndex]->cP(0)
                        - edgePointers[incidentElementsIndex]->cP(1);
    }

    return incidentElementsVectors;
}

bool PatternGeometry::hasAngleSmallerThanThreshold(const std::vector<size_t> &numberOfNodesPerSlot,
                                                   const double &angleThreshold_degrees)
{
    assert(numberOfNodesPerSlot.size() == 7);
    //Initialize helping structs
    if (nodeToSlotMap.empty()) {
        FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlotMap);
    }
    if (correspondingNode.empty()) {
        constructCorrespondingNodeMap(numberOfNodesPerSlot);
    }

    //    const double theta0 = vcg::math::ToDeg(vcg::Point2d(1, 0).Angle());
    //    const double theta1 = vcg::math::ToDeg(vcg::Point2d(1, 1).Angle());
    //    const double theta2 = vcg::math::ToDeg(vcg::Point2d(0, 1).Angle());
    //    const double theta3 = vcg::math::ToDeg(vcg::Point2d(-1, 1).Angle());
    //    const double theta4 = vcg::math::ToDeg(vcg::Point2d(-1, 0).Angle());
    //    const double theta5 = vcg::math::ToDeg(vcg::Point2d(-1, -1).Angle() + 2 * M_PI);
    //    const double theta6 = vcg::math::ToDeg(vcg::Point2d(0, -1).Angle() + 2 * M_PI);
    //    const double theta7 = vcg::math::ToDeg(vcg::Point2d(1, -1).Angle() + 2 * M_PI);

    //    std::set<double> test{theta0, theta1, theta2, theta3, theta4, theta5, theta6, theta7};

    bool hasAngleSmallerThan = false;
    //find tiled incident edges for each node
    using EdgeIndex = int;
    using ThetaWithXAxis = double;
    for (size_t vi = 0; vi < VN(); vi++) {
        std::vector<EdgePointer> incidentElementPointers;
        vcg::edge::VEStarVE(&vert[vi], incidentElementPointers);
        const size_t numberOfIncidentEdges = incidentElementPointers.size();
        if (numberOfIncidentEdges == 0) {
            continue;
        }
        std::vector<VectorType> incidentEdgeVectors
            = getEdgeVectorsWithVertexAsOrigin(incidentElementPointers, vi);

        std::vector<VectorType> tiledIncidentVectors = incidentEdgeVectors;
        const size_t slotIndex = nodeToSlotMap[vi];
        if (slotIndex == 2) {
            //NOTE:I assume that base triangle slots 1,2(bottom triangle nodes) are not used
            std::cerr
                << "Slot of bottom base triangle nodes detected.This case is not handled.Exiting.."
                << std::endl;
            std::terminate();
        } else if (slotIndex == 3 || slotIndex == 5) {
            //NOTE: I don't need to check both triangle edges
            std::vector<PatternGeometry::EdgePointer> correspondingVertexIncidentElementPointers;
            vcg::edge::VEStarVE(&vert[correspondingNode[vi]],
                                correspondingVertexIncidentElementPointers);
            std::vector<VectorType> correspondingVertexIncidentVectors
                = getEdgeVectorsWithVertexAsOrigin(correspondingVertexIncidentElementPointers,
                                                   correspondingNode[vi]);
            //            const CoordType &correspondingVertexPosition = vert[correspondingNode[vi]].cP();
            vcg::Matrix33d R;
            if (slotIndex == 3) {
                R = vcg::RotationMatrix(vcg::Point3d(0, 0, 1), vcg::math::ToRad(-360.0 / fanSize));
            } else {
                R = vcg::RotationMatrix(vcg::Point3d(0, 0, 1), vcg::math::ToRad(360.0 / fanSize));
            }
            //            const CoordType &correspondingVertexPosition_rotated = vert[vi].cP();
            std::transform(
                correspondingVertexIncidentVectors.begin(),
                correspondingVertexIncidentVectors.end(),
                correspondingVertexIncidentVectors.begin(),
                [&](const VectorType &v) {
                    //                               CoordType rotatedTarget = v * R;
                    //                               v = rotatedTarget - correspondingVertexPosition_rotated;
                    return R * v;
                });
            tiledIncidentVectors.insert(tiledIncidentVectors.end(),
                                        correspondingVertexIncidentVectors.begin(),
                                        correspondingVertexIncidentVectors.end());
        } else if (slotIndex == 4) {
            std::vector<VectorType> reversedIncidentElementVectors(incidentEdgeVectors.size());
            std::transform(incidentEdgeVectors.begin(),
                           incidentEdgeVectors.end(),
                           reversedIncidentElementVectors.begin(),
                           [&](const VectorType &v) { return -v; });
            //NOTE: for checking the angles not all opposite incident element vectors are needed but rather the two "extreme" ones of slot 4
            //here I simply add them all
            tiledIncidentVectors.insert(tiledIncidentVectors.end(),
                                        reversedIncidentElementVectors.begin(),
                                        reversedIncidentElementVectors.end());
        }
        //        std::vector<std::array<double, 3>> edgePoints;
        //        for (int tiledVectorIndex = 0; tiledVectorIndex < tiledIncidentVectors.size();
        //             tiledVectorIndex++) {
        //            edgePoints.push_back(
        //                std::array<double, 3>{vert[vi].cP()[0], vert[vi].cP()[1], vert[vi].cP()[2]});
        //            edgePoints.push_back(
        //                std::array<double, 3>{(vert[vi].cP() + tiledIncidentVectors[tiledVectorIndex])[0],
        //                                      (vert[vi].cP() + tiledIncidentVectors[tiledVectorIndex])[1],
        //                                      (vert[vi].cP() + tiledIncidentVectors[tiledVectorIndex])[2]});
        //        }
        //        polyscope::registerCurveNetworkLine("temp", edgePoints);
        //        polyscope::show();
        if (tiledIncidentVectors.size() == 1) {
            continue;
        }

        std::vector<double> thetaAnglesOfIncidentVectors(tiledIncidentVectors.size());
        std::transform(tiledIncidentVectors.begin(),
                       tiledIncidentVectors.end(),
                       thetaAnglesOfIncidentVectors.begin(),
                       [](const VectorType &v) { return vcg::Point2d(v[0], v[1]).Angle(); });
        //sort them using theta angles
        std::sort(thetaAnglesOfIncidentVectors.begin(), thetaAnglesOfIncidentVectors.end());

        //find nodes that contain incident edges with relative angles less than the threshold
        const double angleThreshold_rad = vcg::math::ToRad(angleThreshold_degrees);
        for (int thetaAngleIndex = 1; thetaAngleIndex < thetaAnglesOfIncidentVectors.size();
             thetaAngleIndex++) {
            const double absAngleDifference = std::abs(
                thetaAnglesOfIncidentVectors[thetaAngleIndex]
                - thetaAnglesOfIncidentVectors[thetaAngleIndex - 1]);
            if (absAngleDifference < angleThreshold_rad
                /*&& absAngleDifference > vcg::math::ToRad(0.01)*/) {
                //                std::cout << "Found angDiff:" << absAngleDifference << std::endl;
                //                vert[vi].C() = vcg::Color4b::Magenta;
                //                hasAngleSmallerThan = true;
                return true;
            }
        }
        const double firstLastPairAngleDiff = std::abs(
            thetaAnglesOfIncidentVectors[0]
            - thetaAnglesOfIncidentVectors[thetaAnglesOfIncidentVectors.size() - 1]);
        if (firstLastPairAngleDiff < angleThreshold_rad && firstLastPairAngleDiff > 0.01) {
            return true;
        }
    }
    return false;
}

bool PatternGeometry::hasValenceGreaterThan(const std::vector<size_t> &numberOfNodesPerSlot,
                                            const size_t &valenceThreshold)
{
    //Initialize helping structs
    if (nodeToSlotMap.empty()) {
        FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlotMap);
    }
    if (correspondingNode.empty()) {
        constructCorrespondingNodeMap(numberOfNodesPerSlot);
    }
    if (nodeTiledValence.empty()) {
        constructNodeToTiledValenceMap(numberOfNodesPerSlot);
    }

    //Check tiled valence of the pattern's nodes
    bool hasNodeWithValenceGreaterThanDesired = false;
    for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) {
        const size_t tiledValence = nodeTiledValence[nodeIndex];
        if (tiledValence > valenceThreshold) {
            vert[nodeIndex].C() = vcg::Color4b::LightRed;
            hasNodeWithValenceGreaterThanDesired = true;
        }
    }
    return hasNodeWithValenceGreaterThanDesired;
}

bool PatternGeometry::hasDanglingEdges(
    const std::vector<size_t> &numberOfNodesPerSlot) {
    //Initialize helping structs
    if (nodeToSlotMap.empty()) {
        FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlotMap);
    }
    if (correspondingNode.empty()) {
        constructCorrespondingNodeMap(numberOfNodesPerSlot);
    }
  if (nodeTiledValence.empty()) {
    constructNodeToTiledValenceMap(numberOfNodesPerSlot);
  }

  //Check tiled valence of the pattern's nodes
  bool hasDanglingEdges = false;
  for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) {
    const size_t tiledValence = nodeTiledValence[nodeIndex];
    if (tiledValence == 1) {
      vert[nodeIndex].C() = vcg::Color4b::Red;
      hasDanglingEdges = true;
    }
  }
  return hasDanglingEdges;
}

bool PatternGeometry::hasUntiledDanglingEdges() {
  //  vcg::tri::Clean<TrianglePatternGeometry>::MergeCloseVertex(*this,
  //  0.0000005);
  //  vcg::tri::Allocator<TrianglePatternGeometry>::CompactEveryVector(*this);
  //  vcg::tri::UpdateTopology<TrianglePatternGeometry>::VertexEdge(*this);
  //  vcg::tri::UpdateTopology<TrianglePatternGeometry>::EdgeEdge(*this);
  bool hasDanglingEdges = false;
  for (size_t vi = 0; vi < vn; vi++) {
    std::vector<PatternGeometry::EdgeType *> connectedEdges;
    vcg::edge::VEStarVE(&vert[vi], connectedEdges);
    const size_t nodeValence = connectedEdges.size();
    if (nodeValence == 1) {
      if (vert[vi].C().operator==(vcg::Color4b(vcg::Color4b::Red))) {

      } else {
        vert[vi].C() = vcg::Color4b::Blue;
      }
      hasDanglingEdges = true;
    }
  }
  return hasDanglingEdges;
}

// TODO: The function expects that the numberOfNodesPerSlot follows a specific
// format and that the vertex container was populated in a particular order.
void PatternGeometry::constructCorrespondingNodeMap(const std::vector<size_t> &numberOfNodesPerSlot)
{
    assert(vn != 0 && !nodeToSlotMap.empty() && correspondingNode.empty()
           && numberOfNodesPerSlot.size() == 7);

    for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) {
        const size_t slotIndex = nodeToSlotMap[nodeIndex];
        if (slotIndex == 1) {
            correspondingNode[nodeIndex] = nodeIndex + 1;
        } else if (slotIndex == 2) {
            correspondingNode[nodeIndex] = nodeIndex - 1;
        } else if (slotIndex == 3) {
            const size_t numberOfNodesBefore = nodeIndex
                                               - std::accumulate(numberOfNodesPerSlot.begin(),
                                                                 numberOfNodesPerSlot.begin() + 3,
                                                                 0);
            correspondingNode[nodeIndex] = std::accumulate(numberOfNodesPerSlot.begin(),
                                                           numberOfNodesPerSlot.begin() + 6,
                                                           0)
                                           - 1 - numberOfNodesBefore;
        } else if (slotIndex == 5) {
            const size_t numberOfNodesAfter = std::accumulate(numberOfNodesPerSlot.begin(),
                                                              numberOfNodesPerSlot.begin() + 6,
                                                              0)
                                              - 1 - nodeIndex;
            correspondingNode[nodeIndex] = numberOfNodesAfter
                                           + std::accumulate(numberOfNodesPerSlot.begin(),
                                                             numberOfNodesPerSlot.begin() + 3,
                                                             0);
        }
    }
}

bool PatternGeometry::isFullyConnectedWhenFanned()
{
    assert(vn != 0 /* && !correspondingNode.empty()*/);
    //  TrianglePatternGeometry copyOfPattern(*this);

    //  // If bottom interface nodes have a valence of zero there definetely more than
    //  // a single cc
    //  bool bottomInterfaceIsConnected = false;
    //  for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) {
    //      if (nodeToSlotMap[nodeIndex] == 1 || nodeToSlotMap[nodeIndex] == 4
    //          || nodeToSlotMap[nodeIndex] == 2) {
    //          std::vector<PatternGeometry::EdgeType *> connectedEdges;
    //          vcg::edge::VEStarVE(&vert[nodeIndex], connectedEdges);
    //          const size_t nodeValence = connectedEdges.size();
    //          if (nodeValence != 0) {
    //              bottomInterfaceIsConnected = true;
    //              break;
    //          }
    //      }
    //  }
    //  if (!bottomInterfaceIsConnected) {
    //    return false;
    //  }

    PatternGeometry fanedPattern = createFan(*this);
    vcg::tri::Clean<PatternGeometry>::MergeCloseVertex(fanedPattern, 0.000000005);
    vcg::tri::Allocator<PatternGeometry>::CompactEveryVector(fanedPattern);
    vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(fanedPattern);
    vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(fanedPattern);
    std::vector<std::pair<int, PatternGeometry::EdgePointer>> eCC;
    vcg::tri::Clean<PatternGeometry>::edgeMeshConnectedComponents(fanedPattern, eCC);

    const bool sideInterfaceIsConnected = 1 == eCC.size();
    if (!sideInterfaceIsConnected) {
        return false;
    }

    //  for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) {
    //    if (nodeTiledValence[nodeIndex] == 0)
    //      continue;
    //    const size_t slotIndex = nodeSlot[nodeIndex];
    //    // connect nodes belonging to first or third slot(nodes on the first
    //    // triangle edge)
    //    //      if (nodeTiledValence[nodeIndex] == 0 && slotIndex == 3) {
    //    //        continue;
    //    //      }
    //    if (slotIndex == 1 || slotIndex == 3) {
    //      assert(correspondingNode.count(nodeIndex) != 0);
    //      const size_t correspondingNodeIndex =
    //      correspondingNode[nodeIndex];
    //      std::vector<TrianglePatternGeometry::EdgeType *> starEdges;
    //      vcg::edge::VEStarVE(&vert[nodeIndex], starEdges);
    //      bool containsEdge = false;
    //      for (TrianglePatternGeometry::EdgeType *e : starEdges) {
    //        assert(vcg::tri::Index(*this, e->V(0)) == nodeIndex ||
    //                vcg::tri::Index(*this, e->V(1)) == nodeIndex);
    //        if (vcg::tri::Index(*this, e->V(0)) == nodeIndex) {
    //          if (vcg::tri::Index(*this, e->V(1)) == correspondingNodeIndex) {
    //            containsEdge = true;
    //            break;
    //          }
    //        } else if (vcg::tri::Index(*this, e->V(1)) == nodeIndex) {
    //          if (vcg::tri::Index(*this, e->V(0)) == correspondingNodeIndex) {
    //            containsEdge = true;
    //            break;
    //          }
    //        }
    //      }
    //      if (!containsEdge) {
    //        vcg::tri::Allocator<TrianglePatternGeometry>::AddEdge(
    //            copyOfPattern, nodeIndex, correspondingNodeIndex);
    //      }
    //    } else if (slotIndex == 2 || slotIndex == 5) {
    //      assert(correspondingNode.count(nodeIndex) != 0);
    //    } else {
    //      assert(correspondingNode.count(nodeIndex) == 0);
    //    }
    //  }

    //  std::vector<std::pair<int, TrianglePatternGeometry::EdgePointer>> eCC;
    //  vcg::tri::Clean<TrianglePatternGeometry>::edgeMeshConnectedComponents(
    //      copyOfPattern, eCC);
    //  size_t numberOfCC_edgeBased = eCC.size();
    //  size_t numberOfCC_vertexBased = numberOfCC_edgeBased;
    //  if (numberOfCC_edgeBased == 1) {
    //    vcg::tri::UpdateTopology<TrianglePatternGeometry>::VertexEdge(
    //        copyOfPattern);
    //    vcg::tri::UpdateTopology<TrianglePatternGeometry>::EdgeEdge(copyOfPattern);
    //    vcg::tri::UpdateFlags<TrianglePatternGeometry>::VertexSetV(copyOfPattern);
    //    vcg::tri::UpdateFlags<TrianglePatternGeometry>::VertexClear(copyOfPattern);
    //    for (size_t ei = 0; ei < copyOfPattern.EN(); ei++) {
    //      copyOfPattern.edge[ei].V(0)->SetV();
    //      copyOfPattern.edge[ei].V(1)->SetV();
    //    }

    //    for (size_t vi = 0; vi < copyOfPattern.VN(); vi++) {
    //      if (!copyOfPattern.vert[vi].IsV()) {
    //        numberOfCC_vertexBased++;
    //      }
    //    }
    //    return numberOfCC_vertexBased;
    //  }

    //  return numberOfCC_edgeBased == 1; // TODO: not good
    return true;
}

bool PatternGeometry::hasIntersectingEdges(
    const std::string &patternBinaryRepresentation,
    const std::unordered_map<size_t, std::unordered_set<size_t>>
        &intersectingEdges) {
  bool containsIntersectingEdges = false;
  for (size_t validEdgeIndex = 0;
       validEdgeIndex < patternBinaryRepresentation.size(); validEdgeIndex++) {
      if (patternBinaryRepresentation[validEdgeIndex] == '1'
          && intersectingEdges.contains(validEdgeIndex)) {
          for (auto edgeIndexIt = intersectingEdges.find(validEdgeIndex)->second.begin();
               edgeIndexIt != intersectingEdges.find(validEdgeIndex)->second.end();
               edgeIndexIt++) {
              if (patternBinaryRepresentation[*edgeIndexIt] == '1') {
                  containsIntersectingEdges = true;
                  //            statistics.numberOfIntersectingEdgesOverAllPatterns++;
                  break; // should be uncommented in order to improve
                         //            performance
              }
          }
          if (containsIntersectingEdges) {
              break; // should be uncommented in order to improve performance
          }
      }
  }

  return containsIntersectingEdges;
}

std::unordered_map<size_t, std::unordered_set<size_t>>
PatternGeometry::getIntersectingEdges(
    size_t &numberOfIntersectingEdgePairs) const {
  std::unordered_map<size_t, std::unordered_set<size_t>> intersectingEdges;
  numberOfIntersectingEdgePairs = 0;
  size_t numberOfEdgePairs;
  if (en == 0 || en == 1) {
    numberOfEdgePairs = 0;
  } else {
    numberOfEdgePairs = binomialCoefficient(en, 2);
  }

  for (size_t edgePairIndex = 0; edgePairIndex < numberOfEdgePairs;
       edgePairIndex++) {
    size_t ei0 =
        en - 2 -
        floor(sqrt(-8 * edgePairIndex + 4 * en * (en - 1) - 7) / 2.0 - 0.5);
    size_t ei1 = edgePairIndex + ei0 + 1 - en * (en - 1) / 2 +
                 (en - ei0) * ((en - ei0) - 1) / 2;
    const vcg::Point2d p0(edge[ei0].cP(0)[0], edge[ei0].cP(0)[1]);
    const float epsilon = 0.000005;
    vcg::Box2d bb0(p0 - vcg::Point2d(epsilon, epsilon),
                   p0 + vcg::Point2d(epsilon, epsilon));
    const vcg::Point2d p1(edge[ei0].cP(1)[0], edge[ei0].cP(1)[1]);
    vcg::Box2d bb1(p1 - vcg::Point2d(epsilon, epsilon),
                   p1 + vcg::Point2d(epsilon, epsilon));
    const vcg::Point2d p2(edge[ei1].cP(0)[0], edge[ei1].cP(0)[1]);
    vcg::Box2d bb2(p2 - vcg::Point2d(epsilon, epsilon),
                   p2 + vcg::Point2d(epsilon, epsilon));
    if (bb2.Collide(bb1) || bb2.Collide(bb0))
      continue;
    const vcg::Point2d p3(edge[ei1].cP(1)[0], edge[ei1].cP(1)[1]);
    vcg::Box2d bb3(p3 - vcg::Point2d(epsilon, epsilon),
                   p3 + vcg::Point2d(epsilon, epsilon));
    if (bb3.Collide(bb1) || bb3.Collide(bb0))
      continue;
    const vcg::Segment2d s0(p0, p1);
    const vcg::Segment2d s1(p2, p3);

    vcg::Point2d intersectionPoint;
    const bool edgesIntersect =
        vcg::SegmentSegmentIntersection(s0, s1, intersectionPoint);
    if (edgesIntersect) {
      numberOfIntersectingEdgePairs++;
      intersectingEdges[ei0].insert(ei1);
      intersectingEdges[ei1].insert(ei0);
    }
  }
  return intersectingEdges;
}

PatternGeometry::PatternGeometry(const std::string &filename,
                                 bool addNormalsIfAbsent) {
  if (!std::filesystem::exists(std::filesystem::path(filename))) {
    assert(false);
    std::cerr << "No flat pattern with name " << filename << std::endl;
    return;
  }
  if (!load(filename)) {
    assert(false);
    std::cerr << "File could not be loaded  " << filename << std::endl;
    return;
  }
  if (addNormalsIfAbsent) {
    addNormals();
  }

  vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
  baseTriangleHeight = computeBaseTriangleHeight();
  baseTriangle = computeBaseTriangle();

  updateEigenEdgeAndVertices();
}

double PatternGeometry::computeBaseTriangleHeight() const
{
    return vcg::Distance(vert[0].cP(), vert[interfaceNodeIndex].cP());
}

void PatternGeometry::deleteDanglingVertices()
{
    vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> pu;
    VCGEdgeMesh::deleteDanglingVertices(pu);
    if (!pu.remap.empty()) {
        interfaceNodeIndex
            = pu.remap[interfaceNodeIndex]; //TODO:Could this be automatically be determined?
    }
}

void PatternGeometry::deleteDanglingVertices(
    vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu)
{
    VCGEdgeMesh::deleteDanglingVertices(pu);
}

vcg::Triangle3<double> PatternGeometry::getBaseTriangle() const
{
    return baseTriangle;
}

void PatternGeometry::addNormals()
{
    bool normalsAreAbsent = vert[0].cN().Norm() < 0.000001;
    if (normalsAreAbsent) {
        for (auto &v : vert) {
            v.N() = CoordType(0, 0, 1);
        }
    }
}

vcg::Triangle3<double> PatternGeometry::computeBaseTriangle() const
{
    const double baseTriangleHeight = computeBaseTriangleHeight();
    double bottomEdgeHalfSize = baseTriangleHeight / std::tan(M_PI / 3);
    CoordType patternCoord0 = vert[0].cP();
    CoordType patternBottomRight = vert[interfaceNodeIndex].cP()
                                   + CoordType(bottomEdgeHalfSize, 0, 0);
    CoordType patternBottomLeft = vert[interfaceNodeIndex].cP()
                                  - CoordType(bottomEdgeHalfSize, 0, 0);
    return vcg::Triangle3<double>(patternCoord0, patternBottomLeft, patternBottomRight);
}

PatternGeometry::PatternGeometry(
    const std::vector<size_t> &numberOfNodesPerSlot,
    const std::vector<vcg::Point2i> &edges) {
  add(numberOfNodesPerSlot, edges);
  addNormals();
  baseTriangleHeight = computeBaseTriangleHeight();
  baseTriangle = computeBaseTriangle();
  updateEigenEdgeAndVertices();
}

//TODO: refactor such that it takes as an input a single pattern and tiles into the desired faces without requiring the each face will contain a single pattern
std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
    std::vector<ConstPatternGeometry> &patterns,
    const std::vector<int> &connectToNeighborsVi,
    const VCGPolyMesh &tileInto,
    const std::vector<int> &perSurfaceFacePatternIndices,
    std::vector<size_t> &tileIntoEdgesToTiledVi,
    std::vector<std::vector<size_t>> &perPatternIndexToTiledPatternEdgeIndex)
{
    perPatternIndexToTiledPatternEdgeIndex.resize(patterns.size());
    std::shared_ptr<PatternGeometry> pTiledPattern(new PatternGeometry);
    std::vector<std::vector<int>> tileIntoEdgeToInterfaceVi(
        tileInto.EN()); //ith element contains all the interface nodes that lay on the ith edge of tileInto
    //Compute the barycentric coords of the verts in the base triangle pattern
    std::vector<std::vector<CoordType>> barycentricCoordinates(patterns.size());
    for (size_t patternIndex = 0; patternIndex < patterns.size(); patternIndex++) {
        const PatternGeometry &pattern = patterns[patternIndex];
        const vcg::Triangle3<double> &baseTriangle = pattern.getBaseTriangle();
        barycentricCoordinates[patternIndex].resize(pattern.VN());
        for (int vi = 0; vi < pattern.VN(); vi++) {
            CoordType barycentricCoords_vi;
            vcg::InterpolationParameters<vcg::Triangle3<double>, double>(baseTriangle,
                                                                         pattern.vert[vi].cP(),
                                                                         barycentricCoords_vi);
            barycentricCoordinates[patternIndex][vi] = barycentricCoords_vi;
        }
    }
    VCGTriMesh tileIntoEdgeMesh;

    assert(vcg::tri::HasFEAdjacency(tileInto));
    assert(vcg::tri::HasFVAdjacency(tileInto));
    for (const VCGPolyMesh::FaceType &f : tileInto.face) {
        const int facePatternIndex = perSurfaceFacePatternIndices[tileInto.getIndex(f)];
        if (facePatternIndex == -1) {
            continue;
        }
        CoordType centerOfFace(0, 0, 0);
        for (size_t vi = 0; vi < f.VN(); vi++) {
            centerOfFace = centerOfFace + f.cP(vi);
        }
        centerOfFace /= f.VN();
        vcg::tri::Allocator<VCGTriMesh>::AddVertex(tileIntoEdgeMesh,
                                                   centerOfFace,
                                                   vcg::Color4b::Yellow);

        std::vector<int> firstInFanConnectToNeighbor_vi(connectToNeighborsVi);
        for (int &vi : firstInFanConnectToNeighbor_vi) {
            vi += pTiledPattern->VN();
        }
        ConstPatternGeometry &pattern = patterns[facePatternIndex];

        for (size_t vi = 0; vi < f.VN(); vi++) {
            auto ep = f.FEp(vi);
            assert(vcg::tri::IsValidPointer(tileInto, ep));
            std::vector<vcg::Point3d> meshTrianglePoints{centerOfFace,
                                                         f.cP(vi),
                                                         vi + 1 == f.VN() ? f.cP(0) : f.cP(vi + 1)};
            //            std::vector<vcg::Point3d> meshTrianglePoints{centerOfFace, ep->cP(0), ep->cP(1)};
            auto fit = vcg::tri::Allocator<VCGTriMesh>::AddFace(tileIntoEdgeMesh,
                                                                meshTrianglePoints[0],
                                                                meshTrianglePoints[1],
                                                                meshTrianglePoints[2]);
            //            CoordType faceNormal = ((meshTrianglePoints[1] - meshTrianglePoints[0])
            //                                    ^ (meshTrianglePoints[2] - meshTrianglePoints[0]))
            //                                       .Normalize();
            //TODO: in planar surfaces I should compute the face of triangle
            CoordType faceNormal = f.cN();

            fit->N() = faceNormal;
            PatternGeometry transformedPattern;
            transformedPattern.copy(pattern);
            //Transform the base triangle nodes to the mesh triangle using barycentric coords
            for (int vi = 0; vi < transformedPattern.VN(); vi++) {
                transformedPattern.vert[vi].P() = CoordType(
                    meshTrianglePoints[0] * barycentricCoordinates[facePatternIndex][vi][0]
                    + meshTrianglePoints[1] * barycentricCoordinates[facePatternIndex][vi][1]
                    + meshTrianglePoints[2] * barycentricCoordinates[facePatternIndex][vi][2]);
            }

            for (VertexType &v : transformedPattern.vert) {
                v.N() = faceNormal;
            }

            vcg::tri::Append<PatternGeometry, PatternGeometry>::Remap remap;
            vcg::tri::Append<PatternGeometry, PatternGeometry>::Mesh(*pTiledPattern,
                                                                     transformedPattern,
                                                                     remap);
            for (size_t ei = 0; ei < pattern.EN(); ei++) {
                perPatternIndexToTiledPatternEdgeIndex[facePatternIndex].push_back(remap.edge[ei]);
            }
            const size_t ei = tileInto.getIndex(ep);
            tileIntoEdgeToInterfaceVi[ei].push_back(remap.vert[pattern.interfaceNodeIndex]);
            //Add edges for connecting the desired vertices
            if (!connectToNeighborsVi.empty()) {
                if (vi + 1 == f.VN()) {
                    for (int connectToNeighborIndex = 0;
                         connectToNeighborIndex < connectToNeighborsVi.size();
                         connectToNeighborIndex++) {
                        auto eIt = vcg::tri::Allocator<PatternGeometry>::AddEdge(
                            *pTiledPattern,
                            firstInFanConnectToNeighbor_vi[connectToNeighborIndex],
                            pTiledPattern->VN() - pattern.VN()
                                + connectToNeighborsVi[connectToNeighborIndex]);
                        perPatternIndexToTiledPatternEdgeIndex[facePatternIndex].push_back(
                            pTiledPattern->getIndex(*eIt));
                    }
                }
                if (vi != 0) {
                    for (int connectToNeighborIndex = 0;
                         connectToNeighborIndex < connectToNeighborsVi.size();
                         connectToNeighborIndex++) {
                        auto eIt = vcg::tri::Allocator<PatternGeometry>::AddEdge(
                            *pTiledPattern,
                            pTiledPattern->VN() - 2 * pattern.VN()
                                + connectToNeighborsVi[connectToNeighborIndex],
                            pTiledPattern->VN() - pattern.VN()
                                + connectToNeighborsVi[connectToNeighborIndex]);
                        perPatternIndexToTiledPatternEdgeIndex[facePatternIndex].push_back(
                            pTiledPattern->getIndex(*eIt));
                    }
                }
            }
            //              tiledPattern.updateEigenEdgeAndVertices();
            //              tiledPattern.registerForDrawing();

            //              polyscope::show();
        }
    }
    vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> pu_vertices;
    vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<EdgePointer> pu_edges;
    pTiledPattern->removeDuplicateVertices(pu_vertices, pu_edges);
    //Update perPattern
    //    for (std::vector<size_t> &tiledPatternEdges : perPatternIndexToTiledPatternEdgeIndex) {
    //        for (size_t &ei : tiledPatternEdges) {
    //            ei = pu_edges.remap[ei];
    //        }

    //        auto end = tiledPatternEdges.end();
    //        for (auto it = tiledPatternEdges.begin(); it != end; ++it) {
    //            end = std::remove(it + 1, end, *it);
    //        }

    //        tiledPatternEdges.erase(end, tiledPatternEdges.end());
    //    }
    const size_t sumOfEdgeIndices = std::accumulate(perPatternIndexToTiledPatternEdgeIndex.begin(),
                                                    perPatternIndexToTiledPatternEdgeIndex.end(),
                                                    0,
                                                    [](const size_t &sum,
                                                       const std::vector<size_t> &v) {
                                                        return sum + v.size();
                                                    });

    assert(pTiledPattern->EN() == sumOfEdgeIndices);

    tileIntoEdgesToTiledVi.clear();
    tileIntoEdgesToTiledVi.resize(tileInto.EN());
    for (int ei = 0; ei < tileInto.EN(); ei++) {
        const std::vector<int> interfaceVis = tileIntoEdgeToInterfaceVi[ei];
        if (interfaceVis.empty()) {
            continue;
        }
        assert(interfaceVis.size() == 1 || interfaceVis.size() == 2);
        assert(
            interfaceVis.size() == 1
            || (interfaceVis.size() == 2
                && (pu_vertices.remap[interfaceVis[0]] == std::numeric_limits<size_t>::max()
                    || pu_vertices.remap[interfaceVis[1]] == std::numeric_limits<size_t>::max())));
        tileIntoEdgesToTiledVi[ei] = pu_vertices.remap[interfaceVis[0]];
    }

    pTiledPattern->deleteDanglingVertices();
    vcg::tri::Allocator<PatternGeometry>::CompactEveryVector(*pTiledPattern);
    pTiledPattern->updateEigenEdgeAndVertices();
    pTiledPattern->save();
    return pTiledPattern;
}

bool PatternGeometry::createHoneycombAtom()
{
    VCGEdgeMesh honeycombQuarter;
    const VCGEdgeMesh::CoordType n(0, 0, 1);
    const double H = 0.2;
    const double height = 1.5 * H;
    const double width = 0.2;
    const double theta = 70;
    const double dy = tan(vcg::math::ToRad(90 - theta)) * width / 2;
    vcg::tri::Allocator<VCGEdgeMesh>::AddVertex(honeycombQuarter,
                                                VCGEdgeMesh::CoordType(0, height / 2, 0),
                                                n);
    vcg::tri::Allocator<VCGEdgeMesh>::AddVertex(honeycombQuarter,
                                                VCGEdgeMesh::CoordType(0, H / 2 - dy, 0),
                                                n);
    vcg::tri::Allocator<VCGEdgeMesh>::AddVertex(honeycombQuarter,
                                                VCGEdgeMesh::CoordType(width / 2, H / 2, 0),
                                                n);
    vcg::tri::Allocator<VCGEdgeMesh>::AddVertex(honeycombQuarter,
                                                VCGEdgeMesh::CoordType(width / 2, 0, 0),
                                                n);
    vcg::tri::Allocator<VCGEdgeMesh>::AddEdge(honeycombQuarter, 0, 1);
    vcg::tri::Allocator<VCGEdgeMesh>::AddEdge(honeycombQuarter, 1, 2);
    vcg::tri::Allocator<VCGEdgeMesh>::AddEdge(honeycombQuarter, 2, 3);

    VCGEdgeMesh honeycombAtom;
    // Top right
    vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::MeshCopy(honeycombAtom, honeycombQuarter);
    // Bottom right
    vcg::Matrix44d rotM;
    rotM.SetRotateDeg(180, vcg::Point3d(1, 0, 0));
    vcg::tri::UpdatePosition<VCGEdgeMesh>::Matrix(honeycombQuarter, rotM);
    vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::Mesh(honeycombAtom, honeycombQuarter);
    // Bottom left
    rotM.SetRotateDeg(180, vcg::Point3d(0, 1, 0));
    vcg::tri::UpdatePosition<VCGEdgeMesh>::Matrix(honeycombQuarter, rotM);
    vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::Mesh(honeycombAtom, honeycombQuarter);
    // Top left
    rotM.SetRotateDeg(180, vcg::Point3d(1, 0, 0));
    vcg::tri::UpdatePosition<VCGEdgeMesh>::Matrix(honeycombQuarter, rotM);
    vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::Mesh(honeycombAtom, honeycombQuarter);

    for (VertexType &v : honeycombAtom.vert) {
        v.P()[2] = 0;
    }

    return true;
}

void PatternGeometry::copy(PatternGeometry &copyFrom)
{
    VCGEdgeMesh::copy(copyFrom);
    baseTriangleHeight = copyFrom.getBaseTriangleHeight();
    baseTriangle = copyFrom.getBaseTriangle();
    interfaceNodeIndex = copyFrom.interfaceNodeIndex;
}

void PatternGeometry::scale(const double &desiredBaseTriangleCentralEdgeSize,
                            const int &interfaceNodeIndex)
{
    const double baseTriangleCentralEdgeSize = computeBaseTriangleHeight();
    const double scaleRatio = desiredBaseTriangleCentralEdgeSize / baseTriangleCentralEdgeSize;
    vcg::tri::UpdatePosition<VCGEdgeMesh>::Scale(*this, scaleRatio);
    baseTriangle = computeBaseTriangle();
    baseTriangleHeight = computeBaseTriangleHeight();
}

void PatternGeometry::createFan(const std::vector<int> &connectToNeighborsVi, const size_t &fanSize)
{
    PatternGeometry rotatedPattern;
    vcg::tri::Append<PatternGeometry, PatternGeometry>::MeshCopy(rotatedPattern, *this);
    for (int rotationCounter = 1; rotationCounter < fanSize; rotationCounter++) {
        vcg::Matrix44d R;
        auto rotationAxis = vcg::Point3d(0, 0, 1);
        R.SetRotateDeg(360.0 / fanSize, rotationAxis);
        vcg::tri::UpdatePosition<PatternGeometry>::Matrix(rotatedPattern, R);
        vcg::tri::Append<PatternGeometry, PatternGeometry>::Mesh(*this, rotatedPattern);
        //Add edges for connecting the desired vertices
        if (!connectToNeighborsVi.empty()) {
            if (rotationCounter == fanSize - 1) {
                for (int connectToNeighborIndex = 0;
                     connectToNeighborIndex < connectToNeighborsVi.size();
                     connectToNeighborIndex++) {
                    vcg::tri::Allocator<PatternGeometry>::AddEdge(
                        *this,
                        connectToNeighborsVi[connectToNeighborIndex],
                        this->VN() - rotatedPattern.VN()
                            + connectToNeighborsVi[connectToNeighborIndex]);
                }
            }
            for (int connectToNeighborIndex = 0;
                 connectToNeighborIndex < connectToNeighborsVi.size();
                 connectToNeighborIndex++) {
                vcg::tri::Allocator<PatternGeometry>::AddEdge(
                    *this,
                    this->VN() - 2 * rotatedPattern.VN()
                        + connectToNeighborsVi[connectToNeighborIndex],
                    this->VN() - rotatedPattern.VN() + connectToNeighborsVi[connectToNeighborIndex]);
            }
        }
        removeDuplicateVertices();
        updateEigenEdgeAndVertices();
    }
}

double PatternGeometry::getBaseTriangleHeight() const
{
    return baseTriangleHeight;
}

bool PatternGeometry::isInterfaceConnected(const std::unordered_set<VertexIndex> &interfaceNodes)
{
    std::unordered_set<VertexIndex> interfaceNodesConnected;

    for (int ei = 0; ei < EN(); ei++) {
        const int evi0 = getIndex(edge[ei].cV(0));
        const int evi1 = getIndex(edge[ei].cV(1));
        if (interfaceNodes.contains(evi0) && !interfaceNodes.contains(evi1)) {
            interfaceNodesConnected.insert(evi0);
        } else if (!interfaceNodes.contains(evi0) && interfaceNodes.contains(evi1)) {
            interfaceNodesConnected.insert(evi1);
        }
    }

    if (interfaceNodesConnected.size() != interfaceNodes.size()) {
        return false;
    }
    return true;
}

std::unordered_set<VertexIndex> PatternGeometry::getInterfaceNodes(
    const std::vector<size_t> &numberOfNodesPerSlot)
{
    if (nodeToSlotMap.empty()) {
        FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlotMap);
    }
    std::unordered_set<VertexIndex> interfaceNodes;
    for (std::pair<VertexIndex, size_t> viSlotPair : nodeToSlotMap) {
        if (viSlotPair.second == 4) {
            interfaceNodes.insert(viSlotPair.first);
        }
    }

    return interfaceNodes;
}