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_STAT_TOPOLOGY_H
#define VCL_ALGORITHMS_MESH_STAT_TOPOLOGY_H
 
#include "../sort.h"
 
#include <vclib/mesh.h>
#include <vclib/space/complex.h>
 
#include <map>
#include <set>
#include <stack>
 
namespace vcl {
 
namespace detail {
 
// function called for the Container Cont of the Mesh
template<typename Cont>
void setReferencedVertices(const auto& mesh, auto& refs, uint& nRefs)
{
    // check if the Cont container of the Mesh has vertex references
    if constexpr (comp::HasVertexReferences<typename Cont::ElementType>) {
        // if there are still some vertices non-referenced
        if (nRefs < mesh.vertexNumber()) {
            constexpr uint ELEM_ID = Cont::ElementType::ELEMENT_ID;
            // for eache element of the Cont container
            for (const auto& el : mesh.template elements<ELEM_ID>()) {
                // for each vertex of the element
                for (uint vi : el.vertexIndices()) {
                    // vertex references shoule never be null
                    assert(vi != UINT_NULL);
                    if (!refs[vi]) {
                        // set the vertex as referenced
                        refs[vi] = true;
                        nRefs++;
                    }
                }
            }
        }
    }
}
 
// function called for each container of the Mesh, using variadic templates
template<typename... Cont>
void setReferencedVertices(
    const auto& mesh,
    auto&       refs,
    uint&       nRefs,
    TypeWrapper<Cont...>)
{
    // call the setReferencedVerticesOnVector function for each container of the
    // mesh
    (setReferencedVertices<Cont>(mesh, refs, nRefs), ...);
}
 
// struct to store the information of the wedge texcoords
template<typename WedgeTexCoordType>
struct WedgeTexCoordsInfo
{
    WedgeTexCoordType texCoord;
    uint              texCoordIndex;
 
    bool operator<(const WedgeTexCoordsInfo& other) const
    {
        if (texCoordIndex == other.texCoordIndex) {
            return texCoord < other.texCoord;
        }
        return texCoordIndex < other.texCoordIndex;
    }
};
 
template<FaceMeshConcept MeshType>
std::vector<bool> nonManifoldVerticesVectorBool(const MeshType& m)
    requires HasPerFaceAdjacentFaces<MeshType>
{
    requirePerFaceAdjacentFaces(m);
 
    using FaceType = MeshType::FaceType;
 
    std::vector<bool> nonManifoldVertices(m.vertexContainerSize(), false);
 
    std::vector<uint> TD(m.vertexContainerSize(), 0);
    std::vector<bool> nonManifoldInc(m.vertexContainerSize(), false);
    // First Loop, count how many faces are incident on a vertex and store it in
    // TD, and flag how many vertices are incident on non manifold edges.
    for (const FaceType& f : m.faces()) {
        for (uint i = 0; i < f.vertexNumber(); ++i) {
            TD[m.index(f.vertex(i))]++;
            if (!isFaceManifoldOnEdge(f, i)) {
                nonManifoldInc[m.index(f.vertex(i))]        = true;
                nonManifoldInc[m.index(f.vertexMod(i + 1))] = true;
            }
        }
    }
 
    std::vector<bool> visited(m.vertexContainerSize(), false);
    for (const FaceType& f : m.faces()) {
        for (uint i = 0; i < f.vertexNumber(); ++i) {
            if (!visited[m.index(f.vertex(i))]) {
                visited[m.index(f.vertex(i))] = true;
                MeshPos pos(&f, i);
                uint    starSize = pos.numberOfAdjacentFacesToV();
                if (starSize != TD[m.index(f.vertex(i))])
                    nonManifoldVertices[m.index(f.vertex(i))] = true;
            }
        }
    }
 
    return nonManifoldVertices;
}
 
template<FaceMeshConcept MeshType>
uint numberEdges(
    const MeshType& m,
    uint&           numBoundaryEdges,
    uint&           numNonManifoldEdges)
{
    std::vector<ConstMeshEdgeUtil<MeshType>> edgeVec =
        fillAndSortMeshEdgeUtilVector(m);
 
    uint numEdges       = 0;
    numBoundaryEdges    = 0;
    numNonManifoldEdges = 0;
 
    size_t f_on_cur_edge = 1;
    for (size_t i = 0; i < edgeVec.size(); ++i) {
        if (((i + 1) == edgeVec.size()) || !(edgeVec[i] == edgeVec[i + 1])) {
            ++numEdges;
            if (f_on_cur_edge == 1)
                ++numBoundaryEdges;
            if (f_on_cur_edge > 2)
                ++numNonManifoldEdges;
            f_on_cur_edge = 1;
        }
        else {
            ++f_on_cur_edge;
        }
    }
    return numEdges;
}
 
inline std::list<uint>                             dummyUintList;
inline std::list<std::list<std::pair<uint, uint>>> dummyListOfLists;
inline std::vector<std::pair<uint, uint>>          dummyVectorOfPairs;
 
} // namespace detail
 
uint countPerFaceVertexReferences(const FaceMeshConcept auto& mesh)
{
    using MeshType = decltype(mesh);
 
    uint nRefs = 0;
 
    if constexpr (TriangleMeshConcept<MeshType>) {
        return mesh.faceNumber() * 3;
    }
    else {
        for (const auto& f : mesh.faces()) {
            nRefs += f.vertexNumber();
        }
    }
 
    return nRefs;
}
 
uint largestFaceSize(const FaceMeshConcept auto& mesh)
{
    using FaceType = RemoveRef<decltype(mesh)>::FaceType;
 
    uint maxFaceSize = 0;
 
    if constexpr (FaceType::VERTEX_NUMBER > 0) {
        return FaceType::VERTEX_NUMBER;
    }
    else {
        for (const auto& f : mesh.faces()) {
            maxFaceSize = std::max(maxFaceSize, f.vertexNumber());
        }
    }
 
    return maxFaceSize;
}
 
uint countTriangulatedTriangles(const FaceMeshConcept auto& mesh)
{
    using MeshType = decltype(mesh);
 
    uint nTris = 0;
 
    for (const auto& f : mesh.faces()) {
        nTris += f.vertexNumber() - 2;
    }
 
    return nTris;
}
 
uint largestPerVertexAdjacentVerticesNumber(const MeshConcept auto& mesh)
{
    requirePerVertexAdjacentVertices(mesh);
 
    uint maxAVN = 0;
 
    for (const auto& v : mesh.vertices()) {
        maxAVN = std::max(maxAVN, v.adjVerticesNumber());
    }
 
    return maxAVN;
}
 
template<uint ELEM_ID>
uint largestPerElementAdjacentFacesNumber(const FaceMeshConcept auto& mesh)
{
    requirePerElementComponent<ELEM_ID, CompId::ADJACENT_FACES>(mesh);
 
    uint maxAFN = 0;
 
    for (const auto& e : mesh.template elements<ELEM_ID>()) {
        maxAFN = std::max(maxAFN, e.adjFacesNumber());
    }
 
    return maxAFN;
}
 
uint largestPerVertexAdjacentFacesNumber(const FaceMeshConcept auto& mesh)
{
    return largestPerElementAdjacentFacesNumber<ElemId::VERTEX>(mesh);
}
 
uint largestPerEdgeAdjacentFacesNumber(const FaceMeshConcept auto& mesh)
    requires EdgeMeshConcept<decltype(mesh)>
{
    return largestPerElementAdjacentFacesNumber<ElemId::EDGE>(mesh);
}
 
template<uint ELEM_ID>
uint largestPerElementAdjacentEdgesNumber(const EdgeMeshConcept auto& mesh)
{
    requirePerElementComponent<ELEM_ID, CompId::ADJACENT_EDGES>(mesh);
 
    uint maxAEN = 0;
 
    for (const auto& e : mesh.template elements<ELEM_ID>()) {
        maxAEN = std::max(maxAEN, e.adjEdgesNumber());
    }
 
    return maxAEN;
}
 
uint largestPerVertexAdjacentEdgesNumber(const EdgeMeshConcept auto& mesh)
{
    return largestPerElementAdjacentEdgesNumber<ElemId::VERTEX>(mesh);
}
 
uint largestPerFaceAdjacentEdgesNumber(const EdgeMeshConcept auto& mesh)
    requires FaceMeshConcept<decltype(mesh)>
{
    return largestPerElementAdjacentEdgesNumber<ElemId::FACE>(mesh);
}
 
uint largestPerEdgeAdjacentEdgesNumber(const EdgeMeshConcept auto& mesh)
{
    return largestPerElementAdjacentEdgesNumber<ElemId::EDGE>(mesh);
}
 
template<FaceMeshConcept MeshType>
uint countVerticesToDuplicateByWedgeTexCoords(
    const MeshType&                     mesh,
    std::vector<std::pair<uint, uint>>& vertWedgeMap =
        detail::dummyVectorOfPairs,
    std::list<uint>& vertsToDuplicate = detail::dummyUintList,
    std::list<std::list<std::pair<uint, uint>>>& facesToReassign =
        detail::dummyListOfLists)
{
    vcl::requirePerFaceMaterialIndex(mesh);
    vcl::requirePerFaceWedgeTexCoords(mesh);
 
    using WedgeTexCoordType  = MeshType::FaceType::WedgeTexCoordType;
    using WedgeTexCoordsInfo = detail::WedgeTexCoordsInfo<WedgeTexCoordType>;
 
    // list of faces that reference a wedge texcoord, with the index of the
    // vertex in the face
    using FaceList = std::list<std::pair<uint, uint>>;
 
    vertWedgeMap.resize(mesh.vertexContainerSize());
    vertsToDuplicate.clear();
    facesToReassign.clear();
 
    uint count = 0;
 
    // for each vertex, I'll store a map of WedgeTexCoordsInfo
    // each element of the map represent a unique wedge texcoord(the texcoord
    // itself and the index of the texcoord), and for each element it maps a
    // list of faces that reference the texcoord
    std::vector<std::map<WedgeTexCoordsInfo, FaceList>> wedges(
        mesh.vertexContainerSize());
 
    for (const auto& f : mesh.faces()) {
        for (uint i = 0; i < f.vertexNumber(); ++i) {
            uint vi = f.vertexIndex(i);
 
            // check if the i-th wedge texcoord of the face already exists
            WedgeTexCoordsInfo wi = {f.wedgeTexCoord(i), f.materialIndex()};
            auto               it = wedges[vi].find(wi);
            if (it == wedges[vi].end()) { // if it doesn't exist, add it
                // if there was already a texcoord for the vertex, it means that
                // the vertex will be duplicated
                if (!wedges[vi].empty()) {
                    count++;
                }
                wedges[vi][wi].emplace_back(f.index(), i);
            }
            else {
                // if it exists, add the face to the list of faces that
                // reference the texcoord
                it->second.emplace_back(f.index(), i);
            }
        }
    }
 
    // for each vertex, check if there are multiple texcoords
    // note: here we will modify the maps for convenience: at the end of this
    // loop, the info will be contained in the two lists vertsToDuplicate and
    // facesToReassign (the wedges vector of map will be inconsistent)
    for (uint vi = 0; auto& map : wedges) {
        if (map.size() > 1) {
            // there are multiple texcoords for the vertex vi, so it will be
            // duplicated
 
            // remove from the map the element having the higher number of
            // faces referencing the texcoord (we do this to reassign the less
            // number of faces)
            auto it = std::max_element(
                map.begin(), map.end(), [](const auto& a, const auto& b) {
                    return a.second.size() < b.second.size();
                });
            // store the reference of the texcoord to keep (pair face/vertex
            // index in the face)
            vertWedgeMap[vi] = it->second.front();
            map.erase(it);
 
            // store the vertex to duplicate, and the faces to reassign
            for (auto& [wi, fl] : map) {
                vertsToDuplicate.push_back(vi);
                facesToReassign.push_back(std::move(fl));
            }
        }
        else {
            if (map.size() == 1) {
                // there is only one texcoord for the vertex vi, so store its
                // reference (pair face/vertex index in the face)
                vertWedgeMap[vi] = map.begin()->second.front();
            }
            else { // the vertex was unreferenced
                vertWedgeMap[vi] = {UINT_NULL, UINT_NULL};
            }
        }
        ++vi;
    }
 
    assert(vertWedgeMap.size() == mesh.vertexContainerSize());
    assert(vertsToDuplicate.size() == count);
    assert(facesToReassign.size() == count);
 
    return count;
}
 
template<typename Container, MeshConcept MeshType>
Container referencedVertices(
    const MeshType& mesh,
    uint&           nUnref,
    bool            onlyFaces = false)
{
    using VertexType = MeshType::VertexType;
 
    uint nRefs = 0;
 
    Container refVertices(mesh.vertexContainerSize(), false);
 
    if (onlyFaces) {
        if constexpr (FaceMeshConcept<MeshType>) {
            using FaceContainer = MeshType::FaceContainer;
 
            detail::setReferencedVertices(
                mesh, refVertices, nRefs, TypeWrapper<FaceContainer>());
        }
    }
    else {
        detail::setReferencedVertices(
            mesh, refVertices, nRefs, typename MeshType::Containers());
    }
 
    nUnref = mesh.vertexNumber() - nRefs;
 
    return refVertices;
}
 
template<MeshConcept MeshType>
uint numberUnreferencedVertices(const MeshType& m, bool onlyFaces = false)
{
    uint nV = 0;
    // store the number of unref vertices into nV
    referencedVertices<std::vector<bool>>(m, nV, onlyFaces);
 
    return nV;
}
 
template<FaceMeshConcept MeshType>
uint numberNonManifoldVertices(const MeshType& m)
{
    std::vector<bool> nonManifoldVertices =
        detail::nonManifoldVerticesVectorBool(m);
    return std::count(
        nonManifoldVertices.begin(), nonManifoldVertices.end(), true);
}
 
template<FaceMeshConcept MeshType>
bool isWaterTight(const MeshType& m)
{
    uint numEdgeBorder, numNonManifoldEdges;
    detail::numberEdges(m, numEdgeBorder, numNonManifoldEdges);
    return numEdgeBorder == 0 && numNonManifoldEdges == 0;
}
 
template<FaceMeshConcept MeshType>
uint numberHoles(const MeshType& m) requires HasPerFaceAdjacentFaces<MeshType>
{
    requirePerFaceAdjacentFaces(m);
 
    using VertexType = MeshType::VertexType;
    using FaceType   = MeshType::FaceType;
 
    uint loopNum = 0;
 
    // create a vector of bools to keep track of visited faces.
    std::vector<bool> visitedFaces(m.faceContainerSize(), false);
 
    // Traverse the mesh using a depth-first search algorithm to find all the
    // holes.
    for (const FaceType& f : m.faces()) {
        uint e = 0;
        for (const VertexType* v : f.vertices()) {
            if (!visitedFaces[m.index(f)] && f.adjFace(e) == nullptr) {
                MeshPos<FaceType> startPos(&f, e);
                MeshPos<FaceType> curPos = startPos;
                do {
                    curPos.nextEdgeOnBorderAdjacentToV();
                    curPos.flipVertex();
                    visitedFaces[m.index(curPos.face())] = true;
                } while (curPos != startPos);
                ++loopNum;
            }
            ++e;
        }
    }
    return loopNum;
}
 
template<FaceMeshConcept MeshType>
std::vector<std::set<uint>> connectedComponents(const MeshType& m)
    requires HasPerFaceAdjacentFaces<MeshType>
{
    requirePerFaceAdjacentFaces(m);
 
    using FaceType = MeshType::FaceType;
 
    std::vector<std::set<uint>> cc;
 
    // create a vector of bools to keep track of visited faces.
    std::vector<bool> visitedFaces(m.faceContainerSize(), false);
 
    // create a stack to hold the faces that need to be visited during the
    // depth-first search.
    std::stack<const FaceType*> sf;
 
    // traverse the mesh using a depth-first search algorithm to find the
    // connected components.
    for (const FaceType& f : m.faces()) {
        if (!visitedFaces[m.index(f)]) { // first time I see this face
            visitedFaces[m.index(f)] = true;
 
            // new connected component
            cc.emplace_back();
            std::set<uint>& ccf = cc[cc.size() - 1];
            ccf.insert(m.index(f));
 
            // while the stack is empty, visit the adjacent faces of the top
            // face of the stack
            sf.push(&f);
            while (!sf.empty()) {
                const FaceType* fpt = sf.top();
                // remove the top face and add it to the connected component
                sf.pop();
                ccf.insert(m.index(fpt));
 
                // add the adjacent faces of the current visited in the stack
                for (uint j = 0; j < fpt->vertexNumber(); ++j) {
                    const FaceType* adjf = fpt->adjFace(j);
                    // if there is an adj face and it has not been visited
                    if (adjf != nullptr && !visitedFaces[m.index(adjf)]) {
                        sf.push(adjf);
                        visitedFaces[m.index(adjf)] = true;
                    }
                }
            }
        }
    }
    return cc;
}
 
template<FaceMeshConcept MeshType>
uint numberConnectedComponents(const MeshType& m)
{
    return connectedComponents(m).size();
}
 
} // namespace vcl
 
#endif // VCL_ALGORITHMS_MESH_STAT_TOPOLOGY_H