Mesh generalized winding numbers via The Antipodal Method (SIGGRAPH 2026) (#2540)
* initial libigl-friendly implementation of the Antipodal Method for generalized winding numbers. Only contains the triangle-mesh version and Embree intersector support * removed sse intrinsics code from the intersector * split into separate files as per style guide * moved paper/project references to the entry header
This commit is contained in:
Philip Trettner
@@ -0,0 +1,259 @@
|
||||
// This file is part of libigl, a simple c++ geometry processing library.
|
||||
//
|
||||
// Copyright (C) 2026 Philip Trettner <trettner@shapedcode.com>, Cedric Martens <cedric.martens@umontreal.ca>
|
||||
//
|
||||
// This Source Code Form is subject to the terms of the Mozilla Public License
|
||||
// v. 2.0. If a copy of the MPL was not distributed with this file, You can
|
||||
// obtain one at http://mozilla.org/MPL/2.0/.
|
||||
#ifndef IGL_WINDINGNUMBERANTIPODALSCENE_H
|
||||
#define IGL_WINDINGNUMBERANTIPODALSCENE_H
|
||||
|
||||
#include "PI.h"
|
||||
#include "parallel_for.h"
|
||||
|
||||
#include <Eigen/Core>
|
||||
|
||||
#include <cassert>
|
||||
#include <cmath>
|
||||
#include <cstdint>
|
||||
#include <unordered_map>
|
||||
#include <vector>
|
||||
|
||||
namespace igl
|
||||
{
|
||||
/// Precomputed scene for the Antipodal Method generalized winding number
|
||||
///
|
||||
/// The scene stores only the weighted open-boundary edges of the input
|
||||
/// triangle mesh and a fixed antipodal reference direction `x0`. Closed
|
||||
/// (manifold) meshes have an empty boundary and the fractional term is
|
||||
/// exactly zero.
|
||||
///
|
||||
/// Querying the winding number additionally requires an `Intersector` that
|
||||
/// returns the signed ray-mesh crossing count along `(p, x0)` over the
|
||||
/// original mesh; the scene itself is intersector-agnostic. See
|
||||
/// `igl::embree::EmbreeIntersector::signedIntersectionsRay` for an
|
||||
/// optimized concrete implementation.
|
||||
///
|
||||
/// ### Intersector concept
|
||||
/// A type `I` satisfies the concept when it exposes:
|
||||
/// - `using OriginType = ...;` (3D row vector type)
|
||||
/// - `using DirectionType = ...;` (3D row vector type)
|
||||
/// - `int signedIntersectionsRay(
|
||||
/// OriginType origin, DirectionType direction,
|
||||
/// /* defaulted tnear, tfar, mask */) const;`
|
||||
///
|
||||
/// `winding_number` casts query point/direction to the intersector's types
|
||||
/// at the call site, so a `double`-precision scene against a float-only
|
||||
/// intersector works without an adaptor.
|
||||
template <typename Scalar>
|
||||
class WindingNumberAntipodalScene
|
||||
{
|
||||
public:
|
||||
using Point = Eigen::Matrix<Scalar, 1, 3>;
|
||||
using Direction = Eigen::Matrix<Scalar, 1, 3>;
|
||||
|
||||
private:
|
||||
struct WeightedSeg
|
||||
{
|
||||
Point a;
|
||||
Point b;
|
||||
Scalar w;
|
||||
};
|
||||
|
||||
public:
|
||||
/// Build a scene from a 3D triangle mesh.
|
||||
///
|
||||
/// @param[in] V #V by 3 list of vertex positions
|
||||
/// @param[in] F #F by 3 list of triangle indices
|
||||
/// @param[in] x0 unit reference direction (defaults to a fixed non-axis-aligned vector)
|
||||
template <typename DerivedV, typename DerivedF>
|
||||
WindingNumberAntipodalScene(
|
||||
const Eigen::MatrixBase<DerivedV> & V,
|
||||
const Eigen::MatrixBase<DerivedF> & F,
|
||||
const Direction & x0 = default_x0())
|
||||
: m_x0(x0), m_face_count(static_cast<size_t>(F.rows()))
|
||||
{
|
||||
assert(V.cols() == 3 && "WindingNumberAntipodalScene: only 3D vertex positions are supported");
|
||||
assert(F.cols() == 3 && "WindingNumberAntipodalScene: only triangle meshes are supported");
|
||||
build_boundary_segments(V, F, m_boundary);
|
||||
}
|
||||
|
||||
/// Single-point query: full generalized winding number at `p`
|
||||
/// (fractional + signed integer crossings).
|
||||
template <typename Intersector, typename Derivedp>
|
||||
Scalar winding_number(
|
||||
const Intersector & intersector,
|
||||
const Eigen::MatrixBase<Derivedp> & p) const
|
||||
{
|
||||
const Point pp(static_cast<Scalar>(p(0)),
|
||||
static_cast<Scalar>(p(1)),
|
||||
static_cast<Scalar>(p(2)));
|
||||
const Direction x1 = -m_x0;
|
||||
Scalar area = Scalar(0);
|
||||
for (const auto & ws : m_boundary)
|
||||
{
|
||||
const Point v0 = ws.a - pp;
|
||||
const Point v1 = ws.b - pp;
|
||||
area += ws.w * half_solid_angle_unorm(x1, v0, v1);
|
||||
}
|
||||
const Scalar frac = area / (Scalar(2) * Scalar(igl::PI));
|
||||
|
||||
using IO = typename Intersector::OriginType;
|
||||
using ID = typename Intersector::DirectionType;
|
||||
const IO io(static_cast<typename IO::Scalar>(pp(0)),
|
||||
static_cast<typename IO::Scalar>(pp(1)),
|
||||
static_cast<typename IO::Scalar>(pp(2)));
|
||||
const ID id(static_cast<typename ID::Scalar>(m_x0(0)),
|
||||
static_cast<typename ID::Scalar>(m_x0(1)),
|
||||
static_cast<typename ID::Scalar>(m_x0(2)));
|
||||
const int c = intersector.signedIntersectionsRay(io, id);
|
||||
|
||||
return frac + Scalar(c);
|
||||
}
|
||||
|
||||
/// Batch query, parallelized via `igl::parallel_for`.
|
||||
///
|
||||
/// @param[in] intersector Concept-compatible intersector built over the
|
||||
/// same mesh used to construct the scene. Must
|
||||
/// be safe to query concurrently.
|
||||
/// @param[in] O #O by 3 list of query points
|
||||
/// @param[out] W #O by 1 list of winding numbers
|
||||
template <typename Intersector, typename DerivedO, typename DerivedW>
|
||||
void winding_number(
|
||||
const Intersector & intersector,
|
||||
const Eigen::MatrixBase<DerivedO> & O,
|
||||
Eigen::PlainObjectBase<DerivedW> & W) const
|
||||
{
|
||||
W.resize(O.rows(), 1);
|
||||
|
||||
// Adaptive parallel-for threshold.
|
||||
//
|
||||
// The libigl thread pool has a roughly fixed ~1 ms TOTAL overhead per
|
||||
// parallel_for invocation (not per iteration). So we only spawn the
|
||||
// pool when the WHOLE batch is expected to take ≥ 1 ms.
|
||||
//
|
||||
// Per-query work heuristic:
|
||||
// t_q ≈ 50 * B + 100 * sqrt(F) ns
|
||||
// with B = boundary segment count, F = triangle count. Pool wins once
|
||||
// t_q * O > 10^6 ns ⇒ O > 10^6 / t_q
|
||||
// which is exactly parallel_for's `min_parallel` semantics.
|
||||
//
|
||||
// (This is a rough heuristic and should be revisited once parallel_for becomes lower-overhead)
|
||||
const double t_q_ns =
|
||||
50.0 * static_cast<double>(m_boundary.size()) +
|
||||
100.0 * std::sqrt(static_cast<double>(m_face_count));
|
||||
const size_t min_parallel = static_cast<size_t>(
|
||||
std::ceil(1.0e6 / std::max(t_q_ns, 1.0)));
|
||||
|
||||
igl::parallel_for(O.rows(), [&](const int o)
|
||||
{
|
||||
W(o) = winding_number(intersector, O.row(o));
|
||||
}, min_parallel);
|
||||
}
|
||||
|
||||
/// Reference direction `x0` used to evaluate this scene.
|
||||
const Direction & x0() const { return m_x0; }
|
||||
/// Number of weighted boundary edge segments.
|
||||
size_t num_boundary_segments() const { return m_boundary.size(); }
|
||||
/// Number of triangles in the original mesh.
|
||||
size_t num_faces() const { return m_face_count; }
|
||||
|
||||
/// Default reference direction: normalize(1, sqrt(2), sqrt(3)). A fixed
|
||||
/// non-axis-aligned unit vector; any unit vector works per the paper;
|
||||
/// this choice avoids accidental alignment with axis-aligned geometry.
|
||||
static Direction default_x0()
|
||||
{
|
||||
Direction d(Scalar(1),
|
||||
Scalar(std::sqrt(2.0)),
|
||||
Scalar(std::sqrt(3.0)));
|
||||
return d / d.norm();
|
||||
}
|
||||
|
||||
/// Half the signed solid angle subtended by the spherical triangle
|
||||
/// (x1, v0, v1) at the origin, via the unnormalized
|
||||
/// Van Oosterom-Strackee formula. `x1` is expected to be a unit vector
|
||||
/// (the antipodal "south pole" `-x0`); `v0`, `v1` need not be normalized.
|
||||
/// Callers accumulate per-edge contributions and divide by 2π.
|
||||
static Scalar half_solid_angle_unorm(
|
||||
const Direction & x1, const Point & v0, const Point & v1)
|
||||
{
|
||||
const Scalar l0 = v0.norm();
|
||||
const Scalar l1 = v1.norm();
|
||||
const Scalar num = x1.dot(v0.cross(v1));
|
||||
const Scalar denom = l0 * l1
|
||||
+ l1 * x1.dot(v0)
|
||||
+ l0 * x1.dot(v1)
|
||||
+ v0.dot(v1);
|
||||
return std::atan2(num, denom);
|
||||
}
|
||||
|
||||
/// Extract the open boundary as oriented, weighted edge segments.
|
||||
///
|
||||
/// Each undirected edge {min, max} accumulates +1 for every triangle
|
||||
/// that traverses it as min→max and -1 for every traversal max→min.
|
||||
/// Interior edges of an oriented manifold cancel to zero and are
|
||||
/// dropped. Surviving edges are emitted with positive weight; the
|
||||
/// segment direction is flipped when the net count is negative so the
|
||||
/// weight is always > 0. Non-manifold (≥3 incident triangles per edge)
|
||||
/// is handled the same way. The surviving net count becomes the weight.
|
||||
template <typename DerivedV, typename DerivedF>
|
||||
static void build_boundary_segments(
|
||||
const Eigen::MatrixBase<DerivedV> & V,
|
||||
const Eigen::MatrixBase<DerivedF> & F,
|
||||
std::vector<WeightedSeg> & out)
|
||||
{
|
||||
out.clear();
|
||||
std::unordered_map<std::uint64_t, int> counts;
|
||||
counts.reserve(static_cast<size_t>(F.rows()) * 3);
|
||||
|
||||
const auto pack = [](int lo, int hi) -> std::uint64_t
|
||||
{
|
||||
return (static_cast<std::uint64_t>(static_cast<std::uint32_t>(lo)) << 32)
|
||||
| static_cast<std::uint64_t>(static_cast<std::uint32_t>(hi));
|
||||
};
|
||||
|
||||
for (Eigen::Index t = 0; t < F.rows(); ++t)
|
||||
{
|
||||
const int tri[3] = {
|
||||
static_cast<int>(F(t, 0)),
|
||||
static_cast<int>(F(t, 1)),
|
||||
static_cast<int>(F(t, 2))
|
||||
};
|
||||
for (int e = 0; e < 3; ++e)
|
||||
{
|
||||
const int u = tri[e];
|
||||
const int v = tri[(e + 1) % 3];
|
||||
const int lo = u < v ? u : v;
|
||||
const int hi = u < v ? v : u;
|
||||
const std::uint64_t k = pack(lo, hi);
|
||||
counts[k] += (u < v) ? +1 : -1;
|
||||
}
|
||||
}
|
||||
|
||||
out.reserve(counts.size());
|
||||
for (const auto & kv : counts)
|
||||
{
|
||||
const int sum = kv.second;
|
||||
if (sum == 0) continue;
|
||||
const std::uint32_t lo = static_cast<std::uint32_t>(kv.first >> 32);
|
||||
const std::uint32_t hi = static_cast<std::uint32_t>(kv.first & 0xFFFFFFFFu);
|
||||
const Point a(static_cast<Scalar>(V(lo, 0)),
|
||||
static_cast<Scalar>(V(lo, 1)),
|
||||
static_cast<Scalar>(V(lo, 2)));
|
||||
const Point b(static_cast<Scalar>(V(hi, 0)),
|
||||
static_cast<Scalar>(V(hi, 1)),
|
||||
static_cast<Scalar>(V(hi, 2)));
|
||||
WeightedSeg ws;
|
||||
if (sum > 0) { ws.a = a; ws.b = b; ws.w = static_cast<Scalar>(sum); }
|
||||
else { ws.a = b; ws.b = a; ws.w = static_cast<Scalar>(-sum); }
|
||||
out.push_back(ws);
|
||||
}
|
||||
}
|
||||
|
||||
std::vector<WeightedSeg> m_boundary;
|
||||
Direction m_x0;
|
||||
size_t m_face_count = 0;
|
||||
};
|
||||
}
|
||||
|
||||
#endif
|
||||
@@ -386,3 +386,58 @@ igl::embree::EmbreeIntersector
|
||||
ray.hit.instID[0] = RTC_INVALID_GEOMETRY_ID;
|
||||
ray.hit.primID = RTC_INVALID_GEOMETRY_ID;
|
||||
}
|
||||
|
||||
IGL_INLINE int
|
||||
igl::embree::EmbreeIntersector
|
||||
::signedIntersectionsRay(
|
||||
const OriginType & origin,
|
||||
const DirectionType & direction,
|
||||
float tnear,
|
||||
float tfar,
|
||||
int mask) const
|
||||
{
|
||||
struct query_context
|
||||
{
|
||||
RTCRayQueryContext base; // MUST be first for the reinterpret_cast in the filter
|
||||
int sum;
|
||||
};
|
||||
|
||||
query_context q;
|
||||
rtcInitRayQueryContext(&q.base);
|
||||
q.sum = 0;
|
||||
|
||||
RTCRay ray{};
|
||||
ray.org_x = origin[0];
|
||||
ray.org_y = origin[1];
|
||||
ray.org_z = origin[2];
|
||||
ray.dir_x = direction[0];
|
||||
ray.dir_y = direction[1];
|
||||
ray.dir_z = direction[2];
|
||||
ray.tnear = tnear;
|
||||
ray.tfar = tfar;
|
||||
ray.mask = static_cast<unsigned int>(mask);
|
||||
ray.flags = 0;
|
||||
|
||||
RTCOccludedArguments rargs;
|
||||
rtcInitOccludedArguments(&rargs);
|
||||
rargs.flags = (RTCRayQueryFlags)(
|
||||
RTC_RAY_QUERY_FLAG_COHERENT | RTC_RAY_QUERY_FLAG_INVOKE_ARGUMENT_FILTER);
|
||||
rargs.feature_mask = RTC_FEATURE_FLAG_ALL;
|
||||
rargs.context = &q.base;
|
||||
rargs.filter = +[](RTCFilterFunctionNArguments const* fargs)
|
||||
{
|
||||
assert(fargs->N == 1 && fargs->valid[0]);
|
||||
auto const& fray = reinterpret_cast<RTCRay&>(*fargs->ray);
|
||||
auto const& fhit = reinterpret_cast<RTCHit&>(*fargs->hit);
|
||||
auto const d = fray.dir_x * fhit.Ng_x
|
||||
+ fray.dir_y * fhit.Ng_y
|
||||
+ fray.dir_z * fhit.Ng_z;
|
||||
reinterpret_cast<query_context*>(fargs->context)->sum += d > 0.0f ? +1 : -1;
|
||||
// Reject the hit so traversal continues and we visit every crossing.
|
||||
fargs->valid[0] = 0;
|
||||
};
|
||||
rargs.occluded = nullptr;
|
||||
|
||||
rtcOccluded1(scene, &ray, &rargs);
|
||||
return q.sum;
|
||||
}
|
||||
|
||||
@@ -37,6 +37,8 @@ namespace igl
|
||||
public:
|
||||
typedef Eigen::Matrix<float,Eigen::Dynamic,3> PointMatrixType;
|
||||
typedef Eigen::Matrix<int,Eigen::Dynamic,3> FaceMatrixType;
|
||||
typedef Eigen::RowVector3f OriginType;
|
||||
typedef Eigen::RowVector3f DirectionType;
|
||||
public:
|
||||
EmbreeIntersector();
|
||||
private:
|
||||
@@ -150,6 +152,23 @@ namespace igl
|
||||
Hit<float> &hit,
|
||||
int mask = 0xFFFFFFFF) const;
|
||||
|
||||
/// Signed ray-mesh crossing count along `(origin, direction)`. Visits
|
||||
/// every hit (occluded ray + argument filter rejecting each) and
|
||||
/// accumulates `sign(direction · Ng)` per hit.
|
||||
///
|
||||
/// @param[in] origin ray origin
|
||||
/// @param[in] direction ray direction (need not be normalized)
|
||||
/// @param[in] tnear start of ray segment
|
||||
/// @param[in] tfar end of ray segment
|
||||
/// @param[in] mask a 32 bit mask to identify active geometries
|
||||
/// @return signed crossing count
|
||||
int signedIntersectionsRay(
|
||||
const OriginType & origin,
|
||||
const DirectionType& direction,
|
||||
float tnear = 0.0f,
|
||||
float tfar = std::numeric_limits<float>::infinity(),
|
||||
int mask = 0xFFFFFFFF) const;
|
||||
|
||||
private:
|
||||
|
||||
struct Vertex {float x,y,z,a;};
|
||||
|
||||
@@ -0,0 +1,29 @@
|
||||
// This file is part of libigl, a simple c++ geometry processing library.
|
||||
//
|
||||
// Copyright (C) 2026 Philip Trettner <trettner@shapedcode.com>, Cedric Martens <cedric.martens@umontreal.ca>
|
||||
//
|
||||
// This Source Code Form is subject to the terms of the Mozilla Public License
|
||||
// v. 2.0. If a copy of the MPL was not distributed with this file, You can
|
||||
// obtain one at http://mozilla.org/MPL/2.0/.
|
||||
//
|
||||
// Embree-target satellite that emits explicit template instantiations of
|
||||
// `igl::winding_number_antipodal` against `igl::embree::EmbreeIntersector`.
|
||||
// igl_core does not link embree, so these instantiations cannot live in
|
||||
// include/igl/winding_number_antipodal.cpp. In header-only mode this file
|
||||
// is empty after preprocessing.
|
||||
#include "../winding_number_antipodal.h"
|
||||
|
||||
#ifdef IGL_STATIC_LIBRARY
|
||||
// Pull in the function body so we can instantiate it here. In static-lib
|
||||
// mode the header does not transitively include the body.
|
||||
#include "../winding_number_antipodal.cpp"
|
||||
#include "EmbreeIntersector.h"
|
||||
|
||||
// Mirror of the (V, F, O, W) instantiations in include/igl/winding_number.cpp,
|
||||
// adapted to the antipodal API (3D triangle meshes only, with EmbreeIntersector
|
||||
// as the intersector).
|
||||
template void igl::winding_number_antipodal<igl::embree::EmbreeIntersector, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::embree::EmbreeIntersector const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
|
||||
template void igl::winding_number_antipodal<igl::embree::EmbreeIntersector, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::embree::EmbreeIntersector const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
|
||||
template void igl::winding_number_antipodal<igl::embree::EmbreeIntersector, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::embree::EmbreeIntersector const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> >&);
|
||||
template void igl::winding_number_antipodal<igl::embree::EmbreeIntersector, Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<float, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::embree::EmbreeIntersector const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> >&);
|
||||
#endif
|
||||
@@ -0,0 +1,33 @@
|
||||
// This file is part of libigl, a simple c++ geometry processing library.
|
||||
//
|
||||
// Copyright (C) 2026 Philip Trettner <trettner@shapedcode.com>, Cedric Martens <cedric.martens@umontreal.ca>
|
||||
//
|
||||
// This Source Code Form is subject to the terms of the Mozilla Public License
|
||||
// v. 2.0. If a copy of the MPL was not distributed with this file, You can
|
||||
// obtain one at http://mozilla.org/MPL/2.0/.
|
||||
#include "winding_number_antipodal.h"
|
||||
#include "WindingNumberAntipodalScene.h"
|
||||
|
||||
template <
|
||||
typename Intersector,
|
||||
typename DerivedV,
|
||||
typename DerivedF,
|
||||
typename DerivedO,
|
||||
typename DerivedW>
|
||||
IGL_INLINE void igl::winding_number_antipodal(
|
||||
const Eigen::MatrixBase<DerivedV> & V,
|
||||
const Eigen::MatrixBase<DerivedF> & F,
|
||||
const Intersector & intersector,
|
||||
const Eigen::MatrixBase<DerivedO> & O,
|
||||
Eigen::PlainObjectBase<DerivedW> & W)
|
||||
{
|
||||
using Scalar = typename DerivedV::Scalar;
|
||||
igl::WindingNumberAntipodalScene<Scalar> scene(V, F);
|
||||
scene.winding_number(intersector, O, W);
|
||||
}
|
||||
|
||||
// No explicit template instantiations live in this file: the function is
|
||||
// templated on `Intersector`, which is opaque to the core target. Concrete
|
||||
// instantiations against `igl::embree::EmbreeIntersector` live in
|
||||
// include/igl/embree/winding_number_antipodal.cpp (igl_embree target),
|
||||
// because igl_core does not link embree.
|
||||
@@ -0,0 +1,61 @@
|
||||
// This file is part of libigl, a simple c++ geometry processing library.
|
||||
//
|
||||
// Copyright (C) 2026 Philip Trettner <trettner@shapedcode.com>, Cedric Martens <cedric.martens@umontreal.ca>
|
||||
//
|
||||
// This Source Code Form is subject to the terms of the Mozilla Public License
|
||||
// v. 2.0. If a copy of the MPL was not distributed with this file, You can
|
||||
// obtain one at http://mozilla.org/MPL/2.0/.
|
||||
#ifndef IGL_WINDING_NUMBER_ANTIPODAL_H
|
||||
#define IGL_WINDING_NUMBER_ANTIPODAL_H
|
||||
|
||||
#include "igl_inline.h"
|
||||
|
||||
#include <Eigen/Core>
|
||||
|
||||
namespace igl
|
||||
{
|
||||
/// Generalized winding number via the Antipodal Method
|
||||
/// [Martens, Trettner, Bessmeltsev 2026, "The Antipodal Method: Fast,
|
||||
/// Accurate, and Robust 3D Generalized Winding Numbers", SIGGRAPH 2026,
|
||||
/// https://arxiv.org/abs/2605.01536].
|
||||
///
|
||||
/// Project: https://martenscedric.github.io/academic-page/publications/antipodal_wn.html
|
||||
///
|
||||
/// One-shot convenience: builds a `WindingNumberAntipodalScene<Scalar>`
|
||||
/// (with `Scalar` deduced from `DerivedV`) and runs the batch query against
|
||||
/// the supplied `Intersector`. For repeated queries on the same mesh,
|
||||
/// construct the scene once and call its `winding_number(...)` member
|
||||
/// directly.
|
||||
///
|
||||
/// Mirrors the API shape of `igl::winding_number(V, F, O, W)`.
|
||||
///
|
||||
/// @param[in] V #V by 3 list of vertex positions
|
||||
/// @param[in] F #F by 3 list of triangle indices
|
||||
/// @param[in] intersector Concept-compatible intersector built over (V, F).
|
||||
/// `igl::embree::EmbreeIntersector` satisfies the
|
||||
/// concept and works out of the box. Just call
|
||||
/// its `init(V.cast<float>(), F.cast<int>())`.
|
||||
/// @param[in] O #O by 3 list of query points
|
||||
/// @param[out] W #O by 1 list of winding numbers
|
||||
///
|
||||
/// \see WindingNumberAntipodalScene
|
||||
/// \see igl::embree::EmbreeIntersector::signedIntersectionsRay
|
||||
template <
|
||||
typename Intersector,
|
||||
typename DerivedV,
|
||||
typename DerivedF,
|
||||
typename DerivedO,
|
||||
typename DerivedW>
|
||||
IGL_INLINE void winding_number_antipodal(
|
||||
const Eigen::MatrixBase<DerivedV> & V,
|
||||
const Eigen::MatrixBase<DerivedF> & F,
|
||||
const Intersector & intersector,
|
||||
const Eigen::MatrixBase<DerivedO> & O,
|
||||
Eigen::PlainObjectBase<DerivedW> & W);
|
||||
}
|
||||
|
||||
#ifndef IGL_STATIC_LIBRARY
|
||||
# include "winding_number_antipodal.cpp"
|
||||
#endif
|
||||
|
||||
#endif
|
||||
@@ -0,0 +1,176 @@
|
||||
#include <test_common.h>
|
||||
#include <igl/winding_number_antipodal.h>
|
||||
#include <igl/WindingNumberAntipodalScene.h>
|
||||
#include <igl/embree/EmbreeIntersector.h>
|
||||
#include <igl/winding_number.h>
|
||||
#include <igl/read_triangle_mesh.h>
|
||||
|
||||
#include <Eigen/Core>
|
||||
|
||||
#include <cmath>
|
||||
#include <random>
|
||||
|
||||
TEST_CASE("winding_number_antipodal: closed mesh is integer-valued", "[igl/embree]")
|
||||
{
|
||||
Eigen::MatrixXd V;
|
||||
Eigen::MatrixXi F;
|
||||
igl::read_triangle_mesh(test_common::data_path("decimated-knight.obj"), V, F);
|
||||
|
||||
// Build a grid of query points inside the AABB.
|
||||
const Eigen::RowVector3d mn = V.colwise().minCoeff();
|
||||
const Eigen::RowVector3d mx = V.colwise().maxCoeff();
|
||||
const Eigen::RowVector3d c = 0.5 * (mn + mx);
|
||||
Eigen::MatrixXd P(8, 3);
|
||||
for (int i = 0; i < 8; ++i)
|
||||
{
|
||||
P(i,0) = c(0) + ((i & 1) ? 0.1 : -0.1) * (mx(0)-mn(0));
|
||||
P(i,1) = c(1) + ((i & 2) ? 0.1 : -0.1) * (mx(1)-mn(1));
|
||||
P(i,2) = c(2) + ((i & 4) ? 0.1 : -0.1) * (mx(2)-mn(2));
|
||||
}
|
||||
|
||||
Eigen::VectorXd W_ref;
|
||||
igl::winding_number(V, F, P, W_ref);
|
||||
|
||||
igl::embree::EmbreeIntersector e;
|
||||
e.init(V.cast<float>(), F.cast<int>());
|
||||
|
||||
igl::WindingNumberAntipodalScene<double> scene(V, F);
|
||||
Eigen::VectorXd W;
|
||||
scene.winding_number(e, P, W);
|
||||
|
||||
REQUIRE(W.size() == W_ref.size());
|
||||
for (int i = 0; i < W.size(); ++i)
|
||||
{
|
||||
REQUIRE(std::abs(W(i) - W_ref(i)) < 1e-3);
|
||||
}
|
||||
}
|
||||
|
||||
TEST_CASE("winding_number_antipodal: single triangle has fractional contribution", "[igl/embree]")
|
||||
{
|
||||
// Triangle in the z=0 plane, oriented CCW seen from +z.
|
||||
Eigen::MatrixXd V(3, 3);
|
||||
V << 0, 0, 0,
|
||||
1, 0, 0,
|
||||
0, 1, 0;
|
||||
Eigen::MatrixXi F(1, 3);
|
||||
F << 0, 1, 2;
|
||||
|
||||
igl::embree::EmbreeIntersector e;
|
||||
e.init(V.cast<float>(), F.cast<int>());
|
||||
igl::WindingNumberAntipodalScene<double> scene(V, F);
|
||||
|
||||
// Query points well above and below the triangle, near its centroid.
|
||||
const Eigen::RowVector3d c(1.0/3.0, 1.0/3.0, 0.0);
|
||||
Eigen::MatrixXd P(2, 3);
|
||||
P.row(0) = c + Eigen::RowVector3d(0, 0, 0.5);
|
||||
P.row(1) = c + Eigen::RowVector3d(0, 0, -0.5);
|
||||
|
||||
Eigen::VectorXd W;
|
||||
scene.winding_number(e, P, W);
|
||||
|
||||
// For a single oriented triangle, the GWN is roughly ±(solid angle / 4π).
|
||||
// For a point on the normal axis at height h above the centroid of an
|
||||
// equilateral-ish triangle, |W| stays well under 0.5 but is nonzero, and
|
||||
// the two sides have opposite signs.
|
||||
REQUIRE(std::abs(W(0)) > 1e-3);
|
||||
REQUIRE(std::abs(W(1)) > 1e-3);
|
||||
REQUIRE(W(0) * W(1) < 0.0);
|
||||
}
|
||||
|
||||
TEST_CASE("winding_number_antipodal: random non-manifold soup matches reference", "[igl/embree]")
|
||||
{
|
||||
// Fixed seed: deterministic across runs.
|
||||
std::mt19937 rng(0xA1B2C3D4u);
|
||||
std::uniform_real_distribution<double> coord(-1.0, 1.0);
|
||||
std::uniform_int_distribution<int> vidx(0, 9);
|
||||
|
||||
// 10 random vertices.
|
||||
Eigen::MatrixXd V(10, 3);
|
||||
for (int i = 0; i < V.rows(); ++i)
|
||||
for (int d = 0; d < 3; ++d)
|
||||
V(i, d) = coord(rng);
|
||||
|
||||
// 30 random triangles (degenerate ones, repeated indices, are skipped
|
||||
// and resampled so igl::winding_number's solid_angle stays well-defined).
|
||||
Eigen::MatrixXi F(30, 3);
|
||||
for (int t = 0; t < F.rows(); ++t)
|
||||
{
|
||||
int a, b, c;
|
||||
do { a = vidx(rng); b = vidx(rng); c = vidx(rng); }
|
||||
while (a == b || b == c || a == c);
|
||||
F(t, 0) = a; F(t, 1) = b; F(t, 2) = c;
|
||||
}
|
||||
|
||||
// 100 random query points in a slightly larger box.
|
||||
Eigen::MatrixXd P(100, 3);
|
||||
std::uniform_real_distribution<double> qcoord(-1.5, 1.5);
|
||||
for (int i = 0; i < P.rows(); ++i)
|
||||
for (int d = 0; d < 3; ++d)
|
||||
P(i, d) = qcoord(rng);
|
||||
|
||||
Eigen::VectorXd W_ref;
|
||||
igl::winding_number(V, F, P, W_ref);
|
||||
|
||||
igl::embree::EmbreeIntersector e;
|
||||
e.init(V.cast<float>(), F.cast<int>());
|
||||
|
||||
igl::WindingNumberAntipodalScene<double> scene(V, F);
|
||||
Eigen::VectorXd W;
|
||||
scene.winding_number(e, P, W);
|
||||
|
||||
REQUIRE(W.size() == W_ref.size());
|
||||
for (int i = 0; i < W.size(); ++i)
|
||||
{
|
||||
REQUIRE(std::abs(W(i) - W_ref(i)) < 1e-5);
|
||||
}
|
||||
}
|
||||
|
||||
TEST_CASE("winding_number_antipodal: free-function one-shot matches scene", "[igl/embree]")
|
||||
{
|
||||
Eigen::MatrixXd V;
|
||||
Eigen::MatrixXi F;
|
||||
igl::read_triangle_mesh(test_common::data_path("decimated-knight.obj"), V, F);
|
||||
|
||||
const Eigen::RowVector3d c = 0.5 * (V.colwise().minCoeff() + V.colwise().maxCoeff());
|
||||
Eigen::MatrixXd P(4, 3);
|
||||
P.row(0) = c;
|
||||
P.row(1) = c + Eigen::RowVector3d(0.1, 0.0, 0.0);
|
||||
P.row(2) = c + Eigen::RowVector3d(0.0, 0.1, 0.0);
|
||||
P.row(3) = c + Eigen::RowVector3d(1e3, 0.0, 0.0); // far outside
|
||||
|
||||
igl::embree::EmbreeIntersector e;
|
||||
e.init(V.cast<float>(), F.cast<int>());
|
||||
|
||||
Eigen::VectorXd W;
|
||||
igl::winding_number_antipodal(V, F, e, P, W);
|
||||
|
||||
Eigen::VectorXd W_ref;
|
||||
igl::winding_number(V, F, P, W_ref);
|
||||
|
||||
REQUIRE(W.size() == W_ref.size());
|
||||
for (int i = 0; i < W.size(); ++i)
|
||||
REQUIRE(std::abs(W(i) - W_ref(i)) < 1e-3);
|
||||
}
|
||||
|
||||
TEST_CASE("winding_number_antipodal: float precision compiles", "[igl/embree]")
|
||||
{
|
||||
using V3f = Eigen::Matrix<float, Eigen::Dynamic, 3>;
|
||||
using F3i = Eigen::Matrix<int, Eigen::Dynamic, 3>;
|
||||
V3f V(3, 3);
|
||||
V << 0, 0, 0,
|
||||
1, 0, 0,
|
||||
0, 1, 0;
|
||||
F3i F(1, 3);
|
||||
F << 0, 1, 2;
|
||||
|
||||
igl::embree::EmbreeIntersector e;
|
||||
e.init(V, F);
|
||||
igl::WindingNumberAntipodalScene<float> scene(V, F);
|
||||
|
||||
Eigen::Matrix<float, Eigen::Dynamic, 3> P(1, 3);
|
||||
P << 1.0f/3.0f, 1.0f/3.0f, 0.5f;
|
||||
|
||||
Eigen::VectorXf W;
|
||||
scene.winding_number(e, P, W);
|
||||
REQUIRE(std::abs(W(0)) > 1e-3f);
|
||||
}
|
||||
Reference in New Issue
Block a user