face_topology.h

/*****************************************************************************
 * VCLib                                                                     *
 * Visual Computing Library                                                  *
 *                                                                           *
 * Copyright(C) 2021-2025                                                    *
 * Visual Computing Lab                                                      *
 * ISTI - Italian National Research Council                                  *
 *                                                                           *
 * All rights reserved.                                                      *
 *                                                                           *
 * This program is free software; you can redistribute it and/or modify      *
 * it under the terms of the Mozilla Public License Version 2.0 as published *
 * by the Mozilla Foundation; either version 2 of the License, or            *
 * (at your option) any later version.                                       *
 *                                                                           *
 * This program is distributed in the hope that it will be useful,           *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of            *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the              *
 * Mozilla Public License Version 2.0                                        *
 * (https://www.mozilla.org/en-US/MPL/2.0/) for more details.                *
 ****************************************************************************/
 
#ifndef VCL_ALGORITHMS_MESH_FACE_TOPOLOGY_H
#define VCL_ALGORITHMS_MESH_FACE_TOPOLOGY_H
 
#include <vclib/algorithms/core/polygon/ear_cut.h>
#include <vclib/concepts/mesh.h>
#include <vclib/exceptions/mesh.h>
#include <vclib/mesh/iterators/face/edge_adj_face_iterator.h>
#include <vclib/misc/comparators.h>
#include <vclib/space/complex/mesh_pos.h>
#include <vclib/space/core/polygon.h>
#include <vclib/views/mesh.h>
 
#include <set>
 
namespace vcl {
 
template<FaceMeshConcept MeshType, FaceConcept FaceType>
void addTriangleFacesFromPolygon(
    MeshType&                m,
    FaceType&                f,
    const std::vector<uint>& polygon)
{
    using VertexType = MeshType::VertexType;
    using CoordType  = VertexType::CoordType;
 
    // from the ids, create a polygon of coordinates
    std::vector<CoordType> polCoords(polygon.size());
    for (uint i = 0; i < polygon.size(); ++i) {
        if (polygon[i] >= m.vertexContainerSize()) {
            throw BadVertexIndexException(
                "Index " + std::to_string(polygon[i]) +
                " is out of range in Vertex Container.");
        }
        if (m.vertex(polygon[i]).deleted()) {
            throw BadVertexIndexException(
                "Vertex " + std::to_string(polygon[i]) + " is deleted.");
        }
        polCoords[i] = m.vertex(polygon[i]).coord();
    }
 
    // compute earcut of the polygons
    std::vector<uint> tris = earCut(polCoords);
 
    // faux edges management: create a set of unordered edges of the polygon
    // note: we use indices from 0 to polygon.size() because that are the output
    // indices given by the earcut algorithm
    std::set<std::pair<uint, uint>, UnorderedPairComparator<uint>>
        unorderedEdges;
    for (uint i = 0; i < polygon.size(); ++i)
        unorderedEdges.emplace(i, (i + 1) % (uint) polygon.size());
 
    if constexpr (FaceType::VERTEX_NUMBER < 0) {
        f.resizeVertices(3);
    }
 
    // set the first triangle of the loaded polygon
    for (uint i = 0; i < f.vertexNumber(); ++i) {
        f.setVertex(i, polygon[tris[i]]);
    }
 
    if constexpr (face::HasFaceBitFlags<FaceType>) {
        if (unorderedEdges.find(std::make_pair(tris[0], tris[1])) ==
            unorderedEdges.end())
            f.edgeFaux(0) = true;
        if (unorderedEdges.find(std::make_pair(tris[1], tris[2])) ==
            unorderedEdges.end())
            f.edgeFaux(1) = true;
        if (unorderedEdges.find(std::make_pair(tris[2], tris[0])) ==
            unorderedEdges.end())
            f.edgeFaux(2) = true;
    }
 
    // remaining triangles, need to create more faces in the mesh
    for (uint i = 3; i < tris.size(); i += 3) {
        uint ff = m.addFace();
 
        if constexpr (FaceType::VERTEX_NUMBER < 0) {
            m.face(ff).resizeVertices(3);
        }
 
        for (uint j = 0; j < m.face(ff).vertexNumber(); ++j) {
            m.face(ff).setVertex(j, polygon[tris[i + j]]);
        }
 
        if constexpr (face::HasFaceBitFlags<FaceType>) {
            if (unorderedEdges.find(std::make_pair(tris[i], tris[i + 1])) ==
                unorderedEdges.end())
                m.face(ff).edgeFaux(0) = true;
            if (unorderedEdges.find(std::make_pair(tris[i + 1], tris[i + 2])) ==
                unorderedEdges.end())
                m.face(ff).edgeFaux(1) = true;
            if (unorderedEdges.find(std::make_pair(tris[i + 2], tris[i + 0])) ==
                unorderedEdges.end())
                m.face(ff).edgeFaux(2) = true;
        }
    }
}
 
template<FaceMeshConcept MeshType>
uint addTriangleFacesFromPolygon(MeshType& m, const std::vector<uint>& polygon)
{
    uint fid = m.addFace();
    addTriangleFacesFromPolygon(m, m.face(fid), polygon);
    return fid;
}
 
template<FaceConcept FaceType>
bool isFaceManifoldOnEdge(const FaceType& f, uint edge)
    requires comp::HasAdjacentFaces<FaceType>
{
    // Check if the AdjacentFaces component is available for the given face.
    if (!comp::isAdjacentFacesAvailableOn(f)) {
        throw MissingComponentException(
            "Face has no Adjacent Faces component.");
    }
 
    // Check if the edge is a boundary edge.
    if (f.adjFace(edge) == nullptr) {
        return true;
    }
    else { // Check if the edge is shared by exactly two faces.
        return f.adjFace(edge)->indexOfAdjFace(&f) != UINT_NULL;
    }
}
 
template<FaceConcept FaceType>
bool isFaceEdgeOnBorder(const FaceType& f, uint edge)
    requires comp::HasAdjacentFaces<FaceType>
{
    if (!comp::isAdjacentFacesAvailableOn(f)) {
        throw MissingComponentException(
            "Face has no Adjacent Faces component.");
    }
 
    return f->adjFace(edge) == nullptr;
}
 
template<FaceConcept FaceType>
bool checkFlipEdge(const FaceType& f, uint edge)
    requires comp::HasAdjacentFaces<FaceType>
{
    if (!comp::isAdjacentFacesAvailableOn(f)) {
        throw MissingComponentException(
            "Face has no Adjacent Faces component.");
    }
 
    using VertexType = FaceType::VertexType;
 
    if (f.vertexNumber() > 3)
        return false;
 
    if (isFaceEdgeOnBorder(f, edge))
        return false;
 
    const VertexType* v0 = f.vertex(edge);
    const VertexType* v1 = f.vertexMod(edge + 1);
 
    const FaceType* of = f.adjFace(edge);
    uint            oe = of->indexOfAdjFace(&f);
    assert(oe != UINT_NULL);
 
    // check if the vertices of the edge are the same
    // e.g. the mesh has to be well oriented
    if (of->vertex(oe) != v1 || of->vertexMod(oe + 1) != v0)
        return false;
 
    // check if the flipped edge is already present in the mesh
    // f_v2 and of_v2 are the vertices of the new edge
    const VertexType* f_v2  = f.vertexMod(edge + 2);
    const VertexType* of_v2 = of->vertexMod(oe + 2);
 
    MeshPos<FaceType> pos(&f, f_v2);
    MeshPos<FaceType> startPos = pos;
    // loop in the one ring of f_v2
    do {
        pos.nextEdgeAdjacentToV();
        if (pos.adjVertex() == of_v2)
            return false;
    } while (pos != startPos);
 
    return true;
}
 
template<FaceConcept FaceType>
uint edgeAdjacentFacesNumber(const FaceType& f, uint edge)
    requires comp::HasAdjacentFaces<FaceType>
{
    if (!comp::isAdjacentFacesAvailableOn(f)) {
        throw MissingComponentException(
            "Face has no Adjacent Faces component.");
    }
 
    ConstEdgeAdjFaceIterator<FaceType> begin(f, edge), end;
    uint                               cnt = 0;
    for (auto it = begin; it != end; ++it) {
        ++cnt;
    }
 
    return cnt;
}
 
template<FaceConcept FaceType>
uint faceEdgesOnBorderNumber(const FaceType& f)
    requires comp::HasAdjacentFaces<FaceType>
{
    if (!comp::isAdjacentFacesAvailableOn(f)) {
        throw MissingComponentException(
            "Face has no Adjacent Faces component.");
    }
 
    uint cnt = 0;
    for (uint i = 0; i < f.vertexNumber(); ++i)
        if (isFaceEdgeOnBorder(f, i))
            cnt++;
 
    return cnt;
}
 
template<FaceConcept FaceType>
auto faceDihedralAngleOnEdge(const FaceType& f, uint e)
    requires comp::HasAdjacentFaces<FaceType>
{
    if (!comp::isAdjacentFacesAvailableOn(f)) {
        throw MissingComponentException(
            "Face has no Adjacent Faces component.");
    }
 
    /*
     *     v0 ___________ vf1
     *       |\          |
     *       | e\     f1 |
     *       |    \e1    |
     *       |f     \    |
     *       |        \  |
     *       |__________\|
     *    vf0             v1
     */
 
    assert(f.adjFace(e) != nullptr);
    const FaceType& f1 = *(f.adjFace(e));
 
    uint e1 = f1.indexOfAdjFace(&f);
    assert(e1 != UINT_NULL);
 
    const auto& vf0 = *(f.vertexMod((int) e - 1));
    const auto& vf1 = *(f1.vertexMod(e1 - 1));
 
    auto n0 = faceNormal(f);
    auto n1 = faceNormal(f1);
 
    auto off0 = n0 * vf0.coord();
    auto off1 = n1 * vf1.coord();
 
    auto dist01 = off0 - n0 * vf1.coord();
    auto dist10 = off1 - n1 * vf0.coord();
 
    auto sign = std::abs(dist01) > std::abs(dist10) ? dist01 : dist10;
 
    auto angleRad = n0.angle(n1);
    return sign > 0 ? angleRad : -angleRad;
}
 
template<FaceConcept FaceType>
void detachAdjacentFacesOnEdge(FaceType& f, uint edge)
    requires comp::HasAdjacentFaces<FaceType>
{
    if (!comp::isAdjacentFacesAvailableOn(f)) {
        throw MissingComponentException(
            "Face has no Adjacent Faces component.");
    }
 
    FaceType* nextFace = f.adjFace(edge);
 
    // if nextFace == nullptr there is nothing to do
    // the face is already detached on the edge
    if (nextFace != nullptr) {
        FaceType*                          prevFace;
        ConstEdgeAdjFaceIterator<FaceType> begin(f, edge), end;
        for (auto it = begin; it != end; ++it) {
            prevFace = *it;
        }
 
        if (nextFace == prevFace) { // manifold case
            uint en = nextFace->indexOfAdjFace(&f);
            assert(en != UINT_NULL);
            nextFace->setAdjFace(en, nullptr);
        }
        else { // non manifold case
            // the prev face does not have to point to f anymore, but to
            // nextFace
            uint pn = prevFace->indexOfAdjFace(&f);
            assert(pn != UINT_NULL);
            prevFace->setAdjFace(pn, nextFace);
        }
        f.setAdjFace(edge, nullptr);
    }
}
 
template<FaceConcept FaceType>
void detachFace(FaceType& f) requires comp::HasAdjacentFaces<FaceType>
{
    if (!comp::isAdjacentFacesAvailableOn(f)) {
        throw MissingComponentException(
            "Face has no Adjacent Faces component.");
    }
 
    using VertexType = FaceType::VertexType;
 
    for (uint e = 0; e < f.vertexNumber(); ++e) {
        detachAdjacentFacesOnEdge(f, e);
 
        // if the vertices have adjacent faces
        if constexpr (comp::HasAdjacentFaces<VertexType>) {
            if (comp::isAdjacentFacesAvailableOn(f.vertex(e))) {
                VertexType* v    = f.vertex(e);
                uint        vpos = v->indexOfAdjFace(&f);
                if (vpos != UINT_NULL) {   // may happen if vertex adj faces are
                                           // not initialized / updated
                    v->eraseAdjFace(vpos); // the vertex v has not anymore the
                                           // adjacent face f
                }
            }
        }
    }
}
 
template<FaceConcept FaceType>
std::set<const FaceType*> floodFacePatch(
    const FaceType& seed,
    auto&&          faceSelector)
{
    if (!comp::isAdjacentFacesAvailableOn(seed)) {
        throw MissingComponentException(
            "Face has no Adjacent Faces component.");
    }
 
    std::set<const FaceType*> faces;
    // only faces that satisfy the faceSelector will stay on the stack
    std::vector<const FaceType*> stackFaces;
 
    if (!faceSelector(seed))
        return faces;
 
    faces.insert(&seed);
 
    for (const FaceType* adjacent : seed.adjFaces()) {
        // adding neighbor faces (if faceSelector returns true) to the stack
        if (adjacent && faceSelector(*adjacent))
            stackFaces.push_back(adjacent);
    }
 
    // while there aren't other faces on the stack
    while (stackFaces.size() > 0) {
        const FaceType* fi = stackFaces[stackFaces.size() - 1];
        stackFaces.pop_back(); // pop
        faces.insert(fi);
        for (const FaceType* adjacent : fi->adjFaces()) {
            if (adjacent && faceSelector(*adjacent)) {
                if (faces.find(adjacent) == faces.end())
                    stackFaces.push_back(adjacent);
            }
        }
    }
    return faces;
}
 
} // namespace vcl
 
#endif // VCL_ALGORITHMS_MESH_FACE_TOPOLOGY_H