point_sampling.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_POINT_SAMPLING_H
#define VCL_ALGORITHMS_MESH_POINT_SAMPLING_H
 
#include <vclib/algorithms/mesh/shuffle.h>
#include <vclib/algorithms/mesh/stat.h>
 
#include <vclib/mesh.h>
#include <vclib/space/complex.h>
 
namespace vcl {
 
template<MeshConcept MeshType>
auto allVerticesPointSampling(
    const MeshType&    m,
    std::vector<uint>& birthVertices,
    bool               onlySelected = false)
{
    using VertexType = MeshType::VertexType;
    using PointType  = VertexType::PositionType;
 
    PointSampler<PointType> sampler;
 
    // Determine the number of vertices to sample
    uint n = onlySelected ? vertexSelectionNumber(m) : m.vertexNumber();
 
    // Reserve space in the sampler and birthVertices vectors
    sampler.reserve(n);
    birthVertices.reserve(n);
    birthVertices.clear();
 
    // Loop through all the vertices in the mesh and add them to the sampler if
    // they are selected
    for (const VertexType& v : m.vertices()) {
        if (!onlySelected || v.selected()) {
            sampler.add(v);
            birthVertices.push_back(m.index(v));
        }
    }
    return sampler;
}
 
template<MeshConcept MeshType>
auto allVerticesPointSampling(const MeshType& m, bool onlySelected = false)
{
    std::vector<uint> v;
    return allVerticesPointSampling(m, v, onlySelected);
}
 
template<FaceMeshConcept MeshType>
auto allFacesPointSampling(
    const MeshType&    m,
    std::vector<uint>& birthFaces,
    bool               onlySelected = false)
{
    using PointType = MeshType::VertexType::PositionType;
    using FaceType  = MeshType::FaceType;
 
    PointSampler<PointType> sampler;
 
    // Determine the number of vertices to sample
    uint n = onlySelected ? faceSelectionNumber(m) : m.faceNumber();
 
    // Reserve space in the sampler and birthVertices vectors
    sampler.reserve(n);
    birthFaces.reserve(n);
    birthFaces.clear();
 
    for (const FaceType& f : m.faces()) {
        if (!onlySelected || f.selected()) {
            sampler.add(f);
            birthFaces.push_back(m.index(f));
        }
    }
    return sampler;
}
 
template<FaceMeshConcept MeshType>
auto allFacesPointSampling(const MeshType& m, bool onlySelected = false)
{
    std::vector<uint> v;
    return allFacesPointSampling(m, v, onlySelected);
}
 
template<MeshConcept MeshType>
auto vertexUniformPointSampling(
    const MeshType&     m,
    uint                nSamples,
    std::vector<uint>&  birthVertices,
    bool                onlySelected = false,
    std::optional<uint> seed         = std::nullopt)
{
    using PointType = MeshType::VertexType::PositionType;
 
    uint vn = onlySelected ? vertexSelectionNumber(m) : m.vertexNumber();
 
    if (nSamples >= vn) {
        return allVerticesPointSampling(m, birthVertices, onlySelected);
    }
 
    PointSampler<PointType> ps;
 
    // Reserve space in the sampler and birthVertices vectors
    ps.reserve(vn);
    birthVertices.reserve(vn);
    birthVertices.clear();
 
    std::uniform_int_distribution<uint> dist(0, m.vertexContainerSize() - 1);
 
    std::mt19937 gen = randomGenerator(seed);
 
    std::vector<bool> visited(m.vertexContainerSize(), false);
    uint              nVisited = 0;
 
    while (nVisited < nSamples) {
        uint vi = dist(gen);
        if (!m.vertex(vi).deleted() && !visited[vi] &&
            (!onlySelected || m.vertex(vi).selected())) {
            visited[vi] = true;
            nVisited++;
            ps.add(m.vertex(vi));
            birthVertices.push_back(vi);
        }
    }
 
    return ps;
}
 
template<MeshConcept MeshType>
auto vertexUniformPointSampling(
    const MeshType&     m,
    uint                nSamples,
    bool                onlySelected = false,
    std::optional<uint> seed         = std::nullopt)
{
    std::vector<uint> v;
    return vertexUniformPointSampling(m, nSamples, v, onlySelected, seed);
}
 
template<FaceMeshConcept MeshType>
auto faceUniformPointSampling(
    const MeshType&     m,
    uint                nSamples,
    std::vector<uint>&  birthFaces,
    bool                onlySelected = false,
    std::optional<uint> seed         = std::nullopt)
{
    using PointType = MeshType::VertexType::PositionType;
 
    uint fn = onlySelected ? faceSelectionNumber(m) : m.faceNumber();
 
    if (nSamples >= fn) {
        return allFacesPointSampling(m, birthFaces, onlySelected);
    }
 
    PointSampler<PointType> ps;
 
    // Reserve space in the sampler and birthFaces vectors
    ps.reserve(fn);
    birthFaces.reserve(fn);
    birthFaces.clear();
 
    std::uniform_int_distribution<uint> dist(0, m.faceContainerSize() - 1);
 
    std::mt19937 gen = randomGenerator(seed);
 
    std::vector<bool> visited(m.faceContainerSize(), false);
    uint              nVisited = 0;
 
    while (nVisited < nSamples) {
        uint fi = dist(gen);
        if (!m.face(fi).deleted() && !visited[fi] &&
            (!onlySelected || m.face(fi).selected())) {
            visited[fi] = true;
            nVisited++;
            ps.add(m.face(fi));
            birthFaces.push_back(fi);
        }
    }
 
    return ps;
}
 
template<FaceMeshConcept MeshType>
auto faceUniformPointSampling(
    const MeshType&     m,
    uint                nSamples,
    bool                onlySelected = false,
    std::optional<uint> seed         = std::nullopt)
{
    std::vector<uint> v;
    return faceUniformPointSampling(m, nSamples, v, onlySelected, seed);
}
 
template<MeshConcept MeshType, typename ScalarType>
auto vertexWeightedPointSampling(
    const MeshType&                m,
    const std::vector<ScalarType>& weights,
    uint                           nSamples,
    std::vector<uint>&             birthVertices,
    std::optional<uint>            seed = std::nullopt)
{
    using PointType = MeshType::VertexType::PositionType;
 
    if (nSamples >= m.vertexNumber()) {
        return allVerticesPointSampling(m, birthVertices);
    }
 
    PointSampler<PointType> ps;
    ps.reserve(nSamples);
    birthVertices.reserve(nSamples);
    birthVertices.clear();
 
    std::discrete_distribution<> dist(std::begin(weights), std::end(weights));
 
    std::mt19937 gen = randomGenerator(seed);
 
    std::vector<bool> visited(m.vertexContainerSize(), false);
    uint              nVisited = 0;
 
    while (nVisited < nSamples) {
        uint vi = dist(gen);
        if (vi < m.vertexContainerSize() && !m.vertex(vi).deleted() &&
            !visited[vi]) {
            visited[vi] = true;
            nVisited++;
            ps.add(m.vertex(vi));
            birthVertices.push_back(vi);
        }
    }
 
    return ps;
}
 
template<MeshConcept MeshType, typename ScalarType>
auto vertexWeightedPointSampling(
    const MeshType&                m,
    const std::vector<ScalarType>& weights,
    uint                           nSamples,
    std::optional<uint>            seed = std::nullopt)
{
    std::vector<uint> v;
    return vertexWeightedPointSampling(m, weights, nSamples, v, seed);
}
 
template<FaceMeshConcept MeshType, typename ScalarType>
auto faceWeightedPointSampling(
    const MeshType&                m,
    const std::vector<ScalarType>& weights,
    uint                           nSamples,
    std::vector<uint>&             birthFaces,
    std::optional<uint>            seed = std::nullopt)
{
    using PointType = MeshType::VertexType::PositionType;
 
    if (nSamples >= m.faceNumber()) {
        return allFacesPointSampling(m);
    }
 
    PointSampler<PointType> ps;
 
    // Reserve space in the sampler and birthFaces vectors
    ps.reserve(nSamples);
    birthFaces.reserve(nSamples);
    birthFaces.clear();
 
    std::discrete_distribution<> dist(std::begin(weights), std::end(weights));
 
    std::mt19937 gen = randomGenerator(seed);
 
    std::vector<bool> visited(m.faceContainerSize(), false);
    uint              nVisited = 0;
 
    while (nVisited < nSamples) {
        uint fi = dist(gen);
        if (fi < m.faceContainerSize() && !m.face(fi).deleted() &&
            !visited[fi]) {
            visited[fi] = true;
            nVisited++;
            ps.add(m.face(fi));
            birthFaces.push_back(fi);
        }
    }
 
    return ps;
}
 
template<FaceMeshConcept MeshType, typename ScalarType>
auto faceWeightedPointSampling(
    const MeshType&                m,
    const std::vector<ScalarType>& weights,
    uint                           nSamples,
    std::optional<uint>            seed = std::nullopt)
{
    std::vector<uint> v;
    return faceWeightedPointSampling(m, weights, nSamples, v, seed);
}
 
template<MeshConcept MeshType>
auto vertexQualityWeightedPointSampling(
    const MeshType&     m,
    uint                nSamples,
    std::optional<uint> seed = std::nullopt)
{
    requirePerVertexQuality(m);
 
    using VertexType  = MeshType::VertexType;
    using QualityType = VertexType::QualityType;
 
    std::vector<QualityType> weights;
    weights.resize(m.vertexContainerSize(), 0);
    for (const VertexType& v : m.vertices()) {
        weights[m.index(v)] = v.quality();
    }
 
    return vertexWeightedPointSampling(m, weights, nSamples, seed);
}
 
template<FaceMeshConcept MeshType>
auto faceQualityWeightedPointSampling(
    const MeshType&     m,
    uint                nSamples,
    std::optional<uint> seed = std::nullopt)
{
    requirePerFaceQuality(m);
 
    using FaceType    = MeshType::FaceType;
    using QualityType = FaceType::QualityType;
 
    std::vector<QualityType> weights;
    weights.resize(m.faceContainerSize(), 0);
    for (const FaceType& f : m.faces()) {
        weights[m.index(f)] = f.quality();
    }
 
    return faceWeightedPointSampling(m, weights, nSamples, seed);
}
 
template<FaceMeshConcept MeshType>
auto vertexAreaWeightedPointSampling(
    const MeshType&     m,
    uint                nSamples,
    std::optional<uint> seed = std::nullopt)
{
    using VertexType = MeshType::VertexType;
    using ScalarType = VertexType::ScalarType;
    using FaceType   = MeshType::FaceType;
 
    std::vector<ScalarType> weights(m.vertexContainerSize(), 0);
    std::vector<uint>       cnt(m.vertexContainerSize(), 0);
 
    // for each vertex, store in weights the adjacent faces area and their
    // number
    for (const FaceType& f : m.faces()) {
        ScalarType area = faceArea(f);
        for (const VertexType* v : f.vertices()) {
            weights[m.index(v)] += area;
            cnt[m.index(v)]++;
        }
    }
 
    // divide each area sum by the number of adjacent faces
    for (uint i = 0; i < weights.size(); i++) {
        if (cnt[i] > 0)
            weights[i] /= cnt[i];
    }
 
    // use these weights to create a sapler
    return vertexWeightedPointSampling(m, weights, nSamples, seed);
}
 
template<FaceMeshConcept MeshType>
auto faceAreaWeightedPointSampling(
    const MeshType&     m,
    uint                nSamples,
    std::optional<uint> seed = std::nullopt)
{
    using FaceType = MeshType::FaceType;
 
    std::vector<double> weights(m.faceContainerSize());
 
    for (const FaceType& f : m.faces()) {
        weights[m.index(f)] = faceArea(f);
    }
 
    return faceWeightedPointSampling(m, weights, nSamples, seed);
}
 
template<FaceMeshConcept MeshType>
auto montecarloPointSampling(
    const MeshType&     m,
    uint                nSamples,
    std::vector<uint>&  birthFaces,
    std::optional<uint> seed = std::nullopt)
{
    using VertexType = MeshType::VertexType;
    using PointType  = VertexType::PositionType;
    using ScalarType = PointType::ScalarType;
    using FaceType   = MeshType::FaceType;
    using Interval   = std::pair<ScalarType, const FaceType*>;
 
    FirstElementPairComparator<Interval> comparator;
    PointSampler<PointType>              sampler;
 
    // Reserve space in the sampler and birthFaces vectors
    sampler.reserve(nSamples);
    birthFaces.reserve(nSamples);
    birthFaces.clear();
 
    std::uniform_real_distribution<ScalarType> dist(0, 1);
 
    std::mt19937 gen = randomGenerator(seed);
 
    std::vector<Interval> intervals(m.faceNumber());
    uint                  i    = 0;
    ScalarType            area = 0;
    for (const FaceType& f : m.faces()) {
        area += faceArea(f);
        intervals[i] = std::make_pair(area, &f);
        i++;
    }
 
    ScalarType meshArea = intervals.back().first;
    for (uint i = 0; i < nSamples; i++) {
        ScalarType val = meshArea * dist(gen);
        // lower_bound returns the furthermost iterator i in [first, last) such
        // that, for every iterator j in [first, i), *j < value. E.g. An
        // iterator pointing to the first element "not less than" val, or end()
        // if every element is less than val.
        typename std::vector<Interval>::iterator it = std::lower_bound(
            intervals.begin(),
            intervals.end(),
            std::make_pair(val, nullptr),
            comparator);
 
        sampler.add(
            *it->second,
            randomPolygonBarycentricCoordinate<ScalarType>(
                it->second->vertexNumber(), gen));
        birthFaces.push_back(it->second->index());
    }
 
    return sampler;
}
 
template<FaceMeshConcept MeshType>
auto montecarloPointSampling(
    const MeshType&     m,
    uint                nSamples,
    std::optional<uint> seed = std::nullopt)
{
    std::vector<uint> birthFaces;
    return montecarloPointSampling(m, nSamples, birthFaces, seed);
}
 
template<FaceMeshConcept MeshType>
auto stratifiedMontecarloPointSampling(
    const MeshType&     m,
    uint                nSamples,
    std::optional<uint> seed = std::nullopt)
{
    using FaceType   = MeshType::FaceType;
    using PointType  = MeshType::VertexType::PositionType;
    using ScalarType = PointType::ScalarType;
 
    PointSampler<PointType> ps;
 
    std::mt19937 gen = randomGenerator(seed);
 
    double area              = surfaceArea(m);
    double samplePerAreaUnit = nSamples / area;
    // Montecarlo sampling.
    double floatSampleNum = 0.0;
 
    for (const FaceType& f : m.faces()) {
        // compute # samples in the current face (taking into account of the
        // remainders)
        floatSampleNum += faceArea(f) * samplePerAreaUnit;
        int faceSampleNum = (int) floatSampleNum;
        // for every sample p_i in T...
        for (int i = 0; i < faceSampleNum; i++)
            ps.add(
                f,
                randomPolygonBarycentricCoordinate<ScalarType>(
                    f.vertexNumber(), gen));
        floatSampleNum -= (double) faceSampleNum;
    }
 
    return ps;
}
 
template<FaceMeshConcept MeshType>
auto montecarloPoissonPointSampling(
    const MeshType&     m,
    uint                nSamples,
    std::optional<uint> seed = std::nullopt)
{
    using FaceType   = MeshType::FaceType;
    using PointType  = MeshType::VertexType::PositionType;
    using ScalarType = PointType::ScalarType;
 
    PointSampler<PointType> ps;
 
    std::mt19937 gen = randomGenerator(seed);
 
    ScalarType area              = surfaceArea(m);
    ScalarType samplePerAreaUnit = nSamples / area;
 
    for (const FaceType& f : m.faces()) {
        ScalarType areaT  = faceArea(f);
        int faceSampleNum = poissonRandomNumber(areaT * samplePerAreaUnit, gen);
 
        // for every sample p_i in T...
        for (int i = 0; i < faceSampleNum; i++)
            ps.add(
                f,
                randomPolygonBarycentricCoordinate<ScalarType>(
                    f.vertexNumber(), gen));
    }
 
    return ps;
}
 
template<FaceMeshConcept MeshType, typename ScalarType>
auto vertexWeightedMontecarloPointSampling(
    const MeshType&                m,
    const std::vector<ScalarType>& weights,
    uint                           nSamples,
    double                         variance,
    std::optional<uint>            seed = std::nullopt)
{
    using PointType = MeshType::VertexType::PositionType;
    using FaceType  = MeshType::FaceType;
 
    // lambda function to compute radius weighted area of a face
    auto weightedArea = [](const FaceType&                f,
                           const MeshType&                m,
                           const std::vector<ScalarType>& r) -> ScalarType {
        ScalarType averageQ = 0;
        for (uint i = 0; i < f.vertexNumber(); ++i)
            averageQ += r[m.index(f.vertex(i))];
 
        averageQ /= f.vertexNumber();
        return averageQ * averageQ * faceArea(f);
    };
 
    PointSampler<PointType> ps;
 
    std::mt19937 gen = randomGenerator(seed);
 
    std::vector<ScalarType> radius =
        vertexRadiusFromWeights(m, weights, 1.0, variance, true);
 
    ScalarType wArea = 0;
    for (const FaceType& f : m.faces())
        wArea += weightedArea(f, m, radius);
 
    ScalarType samplePerAreaUnit = nSamples / wArea;
    // Montecarlo sampling.
    double floatSampleNum = 0.0;
    for (const FaceType& f : m.faces()) {
        // compute # samples in the current face (taking into account of the
        // remainders)
        floatSampleNum += weightedArea(f, m, radius) * samplePerAreaUnit;
        uint faceSampleNum = (uint) floatSampleNum;
 
        // for every sample p_i in T...
        for (uint i = 0; i < faceSampleNum; i++)
            ps.add(
                f,
                randomPolygonBarycentricCoordinate<ScalarType>(
                    f.vertexNumber(), gen));
 
        floatSampleNum -= (double) faceSampleNum;
    }
 
    return ps;
}
 
template<FaceMeshConcept MeshType>
auto vertexQualityWeightedMontecarloPointSampling(
    const MeshType&     m,
    uint                nSamples,
    double              variance,
    std::optional<uint> seed = std::nullopt)
{
    requirePerVertexQuality(m);
 
    using VertexType  = MeshType::VertexType;
    using QualityType = VertexType::QualityType;
 
    std::vector<QualityType> weights;
    weights.resize(m.vertexContainerSize(), 0);
    for (const VertexType& v : m.vertices()) {
        weights[m.index(v)] = v.quality();
    }
 
    return vertexWeightedMontecarloPointSampling(
        m, weights, nSamples, variance, seed);
}
 
} // namespace vcl
 
#endif // VCL_ALGORITHMS_MESH_POINT_SAMPLING_H