geometry.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_GEOMETRY_H
#define VCL_ALGORITHMS_MESH_STAT_GEOMETRY_H
 
#include "barycenter.h"
 
#include <vclib/algorithms/mesh/face_topology.h>
#include <vclib/space/complex/mesh_inertia.h>
#include <vclib/space/core/matrix.h>
 
namespace vcl {
 
template<FaceMeshConcept MeshType>
double volume(const MeshType& m)
{
    MeshInertia<MeshType> i(m);
    return i.volume();
}
 
template<FaceMeshConcept MeshType>
double surfaceArea(const MeshType& m)
{
    using FaceType = MeshType::FaceType;
 
    double area = 0;
    for (const FaceType& f : m.faces()) {
        area += faceArea(f);
    }
    return area;
}
 
template<FaceMeshConcept MeshType>
double borderLength(const MeshType& m)
{
    using FaceType = MeshType::FaceType;
 
    double l = 0;
    for (const FaceType& f : m.faces()) {
        for (uint i = 0; i < f.vertexNumber(); ++i) {
            if (f.adjFace(i) == nullptr) {
                l += f.vertex(i)->coord().dist(f.vertexMod(i + 1)->coord());
            }
        }
    }
    return l;
}
 
template<MeshConcept MeshType>
auto covarianceMatrixOfPointCloud(const MeshType& m)
{
    using VertexType = MeshType::VertexType;
    using ScalarType = VertexType::CoordType::ScalarType;
 
    auto bar = barycenter(m);
 
    Matrix33<ScalarType> mm;
    mm.setZero();
    // compute covariance matrix
    for (const VertexType& v : m.vertices()) {
        auto e = v.coord() - bar;
        m += e.outerProduct(e);
    }
    return m;
}
 
template<FaceMeshConcept MeshType>
auto covarianceMatrixOfMesh(const MeshType& m)
{
    using VertexType = MeshType::VertexType;
    using FaceType   = MeshType::FaceType;
    using CoordType  = VertexType::CoordType;
    using ScalarType = CoordType::ScalarType;
 
    CoordType            bar = shellBarycenter(m);
    Matrix33<ScalarType> C;
    C.setZero();
    Matrix33<ScalarType> C0;
    C0.setZero();
    C0(0, 0) = C0(1, 1) = 2.0;
    C0(0, 1) = C0(1, 0) = 1.0;
    C0 *= 1 / 24.0;
    // integral of (x,y,0) in the same triangle
    Eigen::Vector3<ScalarType> x;
    x << 1 / 6.0, 1 / 6.0, 0;
    Matrix33<ScalarType> A; // matrix that bring the vertices to (v1-v0,v2-v0,n)
    Matrix33<ScalarType> DC;
 
    for (const FaceType& f : m.faces()) {
        const CoordType& p0 = f.vertex(0)->coord();
        const CoordType& p1 = f.vertex(1)->coord();
        const CoordType& p2 = f.vertex(2)->coord();
        CoordType        n  = (p1 - p0).cross(p2 - p0);
        double           da = n.norm();
        n /= da * da;
 
        CoordType tmpp = p1 - p0;
        for (uint j = 0; j < 3; j++)
            A(j, 0) = tmpp(j);
        tmpp = p2 - p0;
        for (uint j = 0; j < 3; j++)
            A(j, 1) = tmpp(j);
        for (uint j = 0; j < 3; j++)
            A(j, 2) = n(j);
        Eigen::Vector3<ScalarType> delta;
        tmpp = p0 - bar;
        for (uint j = 0; j < 3; j++)
            delta(j) = tmpp(j);
 
        /* DC is calculated as integral of (A*x+delta) * (A*x+delta)^T over the
         * triangle, where delta = v0-bary */
        DC.setZero();
        DC += A * C0 * A.transpose();
        Matrix33<ScalarType> tmp = (A * x) * delta.transpose();
        DC += tmp + tmp.transpose();
        DC += tmp;
        tmp = delta * delta.transpose();
        DC += tmp * 0.5;
        DC *= da; // the determinant of A is also the double area of f
        C += DC;
    }
    return C;
}
 
template<MeshConcept MeshType, typename ScalarType>
std::vector<ScalarType> vertexRadiusFromWeights(
    const MeshType&                m,
    const std::vector<ScalarType>& weights,
    double                         diskRadius,
    double                         radiusVariance,
    bool                           invert = false)
{
    using VertexType = MeshType::VertexType;
 
    std::vector<ScalarType> radius(m.vertexContainerSize());
    auto minmax = std::minmax_element(weights.begin(), weights.end());
 
    float minRad   = diskRadius;
    float maxRad   = diskRadius * radiusVariance;
    float deltaQ   = *minmax.second - *minmax.first;
    float deltaRad = maxRad - minRad;
    for (const VertexType& v : m.vertices()) {
        ScalarType w = weights[m.index(v)];
        radius[m.index(v)] =
            minRad +
            deltaRad *
                ((invert ? *minmax.second - w : w - *minmax.first) / deltaQ);
    }
 
    return radius;
}
 
template<FaceMeshConcept MeshType>
std::vector<std::pair<uint, uint>> creaseFaceEdges(
    const MeshType& m,
    double          angleRadNeg,
    double          angleRadPos,
    bool            alsoBorderEdges = false)
{
    requirePerFaceAdjacentFaces(m);
 
    std::vector<std::pair<uint, uint>> creaseEdges;
 
    for (const auto& f : m.faces()) {
        for (uint i = 0; i < f.vertexNumber(); ++i) {
            if (f.adjFace(i) == nullptr) {
                // border edge
                if (alsoBorderEdges) {
                    creaseEdges.push_back({f.index(), i});
                }
            }
            else {
                // internal edge
                double angle = faceDihedralAngleOnEdge(f, i);
                if (angle < angleRadNeg || angle > angleRadPos) {
                    creaseEdges.push_back({f.index(), i});
                }
            }
        }
    }
    return creaseEdges;
}
 
} // namespace vcl
 
#endif // VCL_ALGORITHMS_MESH_STAT_GEOMETRY_H