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 (parallel_for)
// 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 actual algorithm is unchanged: it reuses igl::WindingNumberAABB (a
// header-only hierarchical accelerator) so it stays as fast as igl::winding_number
// while producing identical results.
 
#include <Eigen/Core>
 
// clang-format off
#include <wmtk/utils/DisableWarnings.hpp>
#include <igl/winding_number.h>      // must precede WindingNumberAABB.h (its guard needs it)
#include <igl/WindingNumberAABB.h>
#include <wmtk/utils/EnableWarnings.hpp>
// clang-format on
 
#include <wmtk/threading/parallel_for.hpp>
 
#include <algorithm>
#include <cmath>
 
namespace wmtk::utils {
 
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 (same as igl::winding_number for triangle meshes).
    igl::WindingNumberAABB<Eigen::Matrix<double, 1, 3>, Eigen::MatrixXd, Eigen::MatrixXi> hier(
        V,
        F);
    hier.grow();
 
    // hier.winding_number(p) is const and used by igl the same way from parallel_for,
    // so concurrent queries against the shared hierarchy are safe.
    threading::parallel_for(
        threading::range(0, static_cast<size_t>(O.rows())),
        [&](const threading::range& r) {
            for (int o = r.begin(); o < r.end(); ++o) {
                W(o) = hier.winding_number(O.row(o));
            }
        },
        std::max(num_threads, 1));
}
 
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;
 
    threading::parallel_for(
        threading::range(0, static_cast<size_t>(O.rows())),
        [&](const threading::range& r) {
            for (size_t o = r.begin(); o < r.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));
            }
        },
        std::max(num_threads, 1));
}
 
} // namespace wmtk::utils