WindingNumber.hpp

#pragma once
 
// A local copy of libigl's winding-number evaluation, but with the per-query
// parallelism driven by wmtk's own threading framework instead of igl::parallel_for.
// This means the winding number honours the requested `num_threads` (like every other
// parallel section in wmtk) rather than always grabbing all hardware cores.
//
// The 3D algorithm is igl's hierarchy (WindingNumberHierarchy, a deterministic copy of
// igl::WindingNumberAABB), so it is as fast as igl::winding_number; unlike igl's, its
// result does not depend on rand() or on the number of threads, to the last bit.
 
#include <Eigen/Core>
 
#include <wmtk/utils/WindingNumberHierarchy.hpp>
 
// clang-format off
#include <wmtk/utils/DisableWarnings.hpp>
#include <igl/remove_unreferenced.h>
#include <wmtk/utils/EnableWarnings.hpp>
// clang-format on
 
#include <wmtk/threading/dynamic_parallel_for.hpp>
 
#include <algorithm>
#include <cassert>
#include <cmath>
#include <vector>
 
namespace wmtk::utils {
 
inline constexpr size_t kWindingNumberChunk = 256;
 
namespace detail {
 
inline double signed_angle_2d(double ax, double ay, double bx, double by, double px, double py)
{
    // Gather vectors to source and destination, and their lengths.
    double o2A[2] = {px - ax, py - ay};
    double o2B[2] = {px - bx, py - by};
    double o2Al = 0;
    double o2Bl = 0;
    for (int i = 0; i < 2; i++) {
        o2Al += o2A[i] * o2A[i];
        o2Bl += o2B[i] * o2B[i];
    }
    o2Al = std::sqrt(o2Al);
    o2Bl = std::sqrt(o2Bl);
    // Normalize, guarding the degenerate case where P coincides with an endpoint.
    for (int i = 0; i < 2; i++) {
        if (o2Al != 0) {
            o2A[i] /= o2Al;
        }
        if (o2Bl != 0) {
            o2B[i] /= o2Bl;
        }
    }
    constexpr double two_pi = 2.0 * 3.14159265358979323846;
    return -std::atan2(o2B[0] * o2A[1] - o2B[1] * o2A[0], o2B[0] * o2A[0] + o2B[1] * o2A[1]) /
           two_pi;
}
 
} // namespace detail
 
inline double
winding_number_2d_point(const Eigen::MatrixXd& V, const Eigen::MatrixXi& E, double px, double py)
{
    double w = 0;
    for (Eigen::Index f = 0; f < E.rows(); ++f) {
        const Eigen::Index a = E(f, 0);
        const Eigen::Index b = E(f, 1);
        w += detail::signed_angle_2d(V(a, 0), V(a, 1), V(b, 0), V(b, 1), px, py);
    }
    return w;
}
 
inline void winding_number(
    const Eigen::MatrixXd& V,
    const Eigen::MatrixXi& F,
    const Eigen::MatrixXd& O,
    Eigen::VectorXd& W,
    int num_threads)
{
    W.setZero(O.rows());
    if (O.rows() == 0 || F.rows() == 0 || V.rows() == 0) return;
 
    // Build the accelerator once (as igl::winding_number does for triangle meshes).
    const WindingNumberHierarchy hier(V, F);
 
    // hier.winding_number(p) is const and touches no shared state, so concurrent queries
    // against the one hierarchy are safe.
    //
    // Handed out in chunks on demand rather than one fixed slice per thread: a query costs
    // far more near the surface, where the hierarchy is descended to its leaves, than away
    // from it, and the queries of a mesh are not in random order -- the cells of an
    // arrangement cluster by location -- so fixed slices left most threads waiting on the one
    // that drew the surface. Each query is evaluated on its own, so W does not depend on the
    // schedule or the number of threads.
    threading::dynamic_parallel_for(
        static_cast<size_t>(O.rows()),
        std::max(num_threads, 1),
        kWindingNumberChunk,
        [&](const size_t begin, const size_t end) {
            for (size_t o = begin; o < end; ++o) {
                W(o) = hier.winding_number(O.row(o));
            }
        });
}
 
inline bool winding_number_inside(double w)
{
    return w > 0.5;
}
 
inline bool any_winding_number_inside(const Eigen::VectorXd& W)
{
    return (W.array() > 0.5).any();
}
 
inline bool orient_winding_number(Eigen::VectorXd& W)
{
    if (W.size() == 0 || any_winding_number_inside(W)) return false;
    W = -W;
    return true;
}
 
inline void winding_number_by_group(
    const Eigen::MatrixXd& V,
    const Eigen::MatrixXi& F,
    const std::vector<int>& group,
    int n_groups,
    const Eigen::MatrixXd& O,
    Eigen::MatrixXd& W,
    int num_threads)
{
    assert(group.size() == size_t(F.rows()));
    W.setZero(O.rows(), n_groups);
    std::vector<std::vector<int>> faces(n_groups);
    for (Eigen::Index f = 0; f < F.rows(); ++f) {
        assert(group[f] >= 0 && group[f] < n_groups);
        faces[group[f]].push_back(int(f));
    }
    for (int g = 0; g < n_groups; ++g) {
        if (faces[g].empty()) continue;
        Eigen::MatrixXi Fg(faces[g].size(), 3);
        for (size_t i = 0; i < faces[g].size(); ++i) Fg.row(i) = F.row(faces[g][i]);
        // Only the group's own vertices, so that its hierarchy's root box is the group's.
        Eigen::MatrixXd Vg;
        Eigen::MatrixXi Fg_compact;
        Eigen::VectorXi I;
        igl::remove_unreferenced(V, Fg, Vg, Fg_compact, I);
        Eigen::VectorXd w;
        winding_number(Vg, Fg_compact, O, w, num_threads);
        W.col(g) = w;
    }
}
 
inline void winding_number_2d(
    const Eigen::MatrixXd& V,
    const Eigen::MatrixXi& E,
    const Eigen::MatrixXd& O,
    Eigen::VectorXd& W,
    int num_threads)
{
    W.setZero(O.rows());
    if (O.rows() == 0 || E.rows() == 0 || V.rows() == 0) return;
 
    // In chunks on demand, as the 3D winding_number: see there.
    threading::dynamic_parallel_for(
        static_cast<size_t>(O.rows()),
        std::max(num_threads, 1),
        kWindingNumberChunk,
        [&](const size_t begin, const size_t end) {
            for (size_t o = begin; o < end; ++o) {
                const Eigen::Index i = static_cast<Eigen::Index>(o);
                W(i) = winding_number_2d_point(V, E, O(i, 0), O(i, 1));
            }
        });
}
 
} // namespace wmtk::utils