Fast/robust point-in-mesh test: c_pointsinside(method=)

DESCRIPTION

@@ -3,11 +3,14 @@ Package: natcpp
 Title: Fast C++ Primitives for the 'NeuroAnatomy Toolbox'
 Version: 0.3.1.9000
 Authors@R: 
-    person(given = "Gregory",
-           family = "Jefferis",
-           role = c("aut", "cre"),
-           email = "jefferis@gmail.com",
-           comment = c(ORCID = "0000-0002-0587-9355"))
+    c(person(given = "Gregory",
+             family = "Jefferis",
+             role = c("aut", "cre"),
+             email = "jefferis@gmail.com",
+             comment = c(ORCID = "0000-0002-0587-9355")),
+      person(given = "libigl contributors",
+             role = c("ctb", "cph"),
+             comment = "bundled libigl fast winding number code (MPL-2.0); see src/vendor/README.md"))
 Description: Fast functions implemented in C++ via 'Rcpp' to support the
     'NeuroAnatomy Toolbox' ('nat') ecosystem. These functions provide large
     speed-ups for basic manipulation of neuronal skeletons over pure R
@@ -27,6 +30,7 @@ Suggests:
     testthat (>= 3.0.0)
 LinkingTo:
     Rcpp,
+    RcppEigen,
     RcppThread
 Config/testthat/edition: 3
 Encoding: UTF-8

NAMESPACE

@@ -5,6 +5,7 @@ export(c_ListofMatrixRows)
 export(c_coords21dindex)
 export(c_ijkpos)
 export(c_listlengths)
+export(c_pointsinside)
 export(c_seglengths)
 export(c_sub2ind)
 export(c_topntail)

NEWS.md

@@ -5,6 +5,23 @@
   thread policy (respecting `getOption("Ncpus")` and `OMP_THREAD_LIMIT`, else a
   conservative 2) rather than a hard-coded 4. Pass `threads = 0` for all cores,
   or an integer to override. No new package dependency.
+* add `c_pointsinside()`, a robust point-in-mesh test based on the generalised
+  (solid-angle) winding number. Unlike a closest-point signed-distance test it
+  does not depend on surface normals, so it avoids the spurious "outside point
+  classified as inside" results that normal-based tests can give near thin
+  protrusions or sharp features. `threads = NULL` applies the package thread
+  policy.
+* `c_pointsinside(method=)` selects the winding-number back end: `"bruteforce"`
+  is a self-contained O(P*F) test parallelised over points with RcppThread;
+  `"bvh"` uses libigl's "Fast Winding Numbers for Soups and Clouds" (Barill et
+  al. 2018), building a bounding-volume hierarchy once and evaluating each query
+  point in O(log F) so it scales to millions of points on large meshes. The
+  default `"auto"` picks `"bvh"` for large meshes and `"bruteforce"` otherwise;
+  the two agree to within the winding-number tolerance. libigl (MPL-2.0) is
+  vendored under `src/vendor/igl` and requires `RcppEigen`; its bundled Houdini
+  HDK amalgamation is marked a system header (one-line `#pragma`) so its
+  third-party compiler warnings do not surface as install-time `R CMD check`
+  warnings.
 
 # natcpp 0.3.1
 

R/RcppExports.R

@@ -57,6 +57,14 @@ c_coords21dindex <- function(xyz, origin, voxdims, dims, clamp = FALSE) {
     .Call(`_natcpp_c_coords21dindex`, xyz, origin, voxdims, dims, clamp)
 }
 
+c_fast_mesh_winding_number <- function(points, vertices, faces, threads = 4L, accuracy = 2.0) {
+    .Call(`_natcpp_c_fast_mesh_winding_number`, points, vertices, faces, threads, accuracy)
+}
+
+c_mesh_winding_number <- function(points, vertices, faces, threads = 4L) {
+    .Call(`_natcpp_c_mesh_winding_number`, points, vertices, faces, threads)
+}
+
 #' Convert a matrix into list of row vectors
 #'
 #' @details Typically this will be for 3D coordinates but there are no limits

R/inside_mesh.R

@@ -0,0 +1,81 @@
+#' Test which points lie inside a triangle mesh (generalised winding number)
+#'
+#' @description Robust point-in-mesh test based on the generalised (solid-angle)
+#'   winding number. For a closed mesh the winding number is approximately
+#'   \eqn{\pm 1} for interior points and \eqn{0} for exterior points, so
+#'   \code{abs(w) > 0.5} classifies points as inside. Unlike a closest-point
+#'   signed-distance test it does not depend on surface normals and has no
+#'   ray-casting tie-breaking, so it does not produce the spurious "outside
+#'   point classified as inside" results that normal-based tests can give near
+#'   thin protrusions or sharp features.
+#'
+#' @details The mesh should be closed (watertight) and triangular; the result is
+#'   independent of face orientation (winding). It is intended as the
+#'   accelerated back end for \code{nat::pointsinside()}.
+#'
+#'   Two back ends are available, selected by \code{method}:
+#'   \describe{
+#'     \item{\code{"bruteforce"}}{A self-contained \eqn{O(P \times F)}
+#'       implementation (P points, F faces), parallelised over points with
+#'       \pkg{RcppThread}. No setup cost, so it is fastest for small meshes.}
+#'     \item{\code{"bvh"}}{libigl's "Fast Winding Numbers for Soups and Clouds"
+#'       (Barill et al. 2018): a bounding-volume hierarchy is built once over
+#'       the mesh and each query point is then evaluated in \eqn{O(\log F)}, so
+#'       it scales to millions of points on meshes of tens of thousands of
+#'       faces. \code{accuracy} tunes the multipole approximation.}
+#'   }
+#'   \code{"auto"} (the default) picks \code{"bvh"} only for large meshes queried
+#'   by enough points to amortise building the hierarchy, and \code{"bruteforce"}
+#'   otherwise (so a simple mesh such as a cuboid always uses brute force). The
+#'   two back ends agree to within the winding-number tolerance.
+#'
+#' @param points An Nx3 matrix of query point coordinates (or anything
+#'   coercible with \code{as.matrix}).
+#' @param vertices An Nx3 matrix of mesh vertex coordinates.
+#' @param faces An Nx3 integer matrix of 1-based vertex indices (one triangle
+#'   per row), e.g. \code{t(mesh$it)} for an \pkg{rgl} \code{mesh3d}.
+#' @param method Winding-number back end: \code{"auto"} (default), \code{"bvh"}
+#'   or \code{"bruteforce"}. See \strong{Details}.
+#' @param threads Number of threads for parallel computation. The default
+#'   \code{NULL} applies the package thread policy (respecting
+#'   \code{getOption("Ncpus")} and the \code{OMP_THREAD_LIMIT} environment
+#'   variable, else 2). Set to 0 to use all available cores.
+#' @param accuracy libigl accuracy-scale parameter for \code{method = "bvh"}
+#'   (default 2); ignored by the brute-force back end.
+#' @return A logical vector of length \code{nrow(points)} (\code{TRUE} = inside).
+#' @export
+#' @examples
+#' # tetrahedron
+#' V <- rbind(c(0,0,0), c(1,0,0), c(0,1,0), c(0,0,1))
+#' F <- rbind(c(1,3,2), c(1,2,4), c(1,4,3), c(2,3,4))
+#' c_pointsinside(rbind(c(.2,.2,.2), c(2,2,2)), V, F)  # TRUE FALSE
+c_pointsinside <- function(points, vertices, faces,
+                           method = c("auto", "bvh", "bruteforce"),
+                           threads = NULL, accuracy = 2) {
+  method <- match.arg(method)
+  points <- as.matrix(points)
+  vertices <- as.matrix(vertices)
+  faces <- matrix(as.integer(faces), ncol = 3L)
+  threads <- natcpp_threads(threads)
+
+  if (method == "auto") {
+    nf <- nrow(faces)
+    # The libigl BVH has a fixed build cost (~O(F log F)); above ~1000 faces it
+    # amortises quickly and then far outperforms the brute-force O(P*F) test
+    # (measured ~4x at 7k faces, ~25x at 50k). A mistaken "bvh" choice only ever
+    # costs that bounded build, whereas a mistaken "bruteforce" on a large mesh
+    # costs seconds, so we lean towards bvh once the mesh is non-trivial --
+    # except for tiny total workloads, where brute force is instant and needs no
+    # build (a cuboid, F=12, always lands in brute force).
+    method <- if (nf >= 1000L && as.double(nrow(points)) * nf >= 1e6)
+      "bvh" else "bruteforce"
+  }
+
+  w <- if (method == "bvh")
+    c_fast_mesh_winding_number(points, vertices, faces, threads = threads,
+                               accuracy = accuracy)
+  else
+    c_mesh_winding_number(points, vertices, faces, threads = threads)
+
+  abs(w) > 0.5
+}

src/RcppExports.cpp

@@ -1,6 +1,7 @@
 // Generated by using Rcpp::compileAttributes() -> do not edit by hand
 // Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
+#include <RcppEigen.h>
 #include <RcppThread.h>
 #include <Rcpp.h>
 
@@ -81,6 +82,35 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// c_fast_mesh_winding_number
+NumericVector c_fast_mesh_winding_number(NumericMatrix points, NumericMatrix vertices, IntegerMatrix faces, int threads, double accuracy);
+RcppExport SEXP _natcpp_c_fast_mesh_winding_number(SEXP pointsSEXP, SEXP verticesSEXP, SEXP facesSEXP, SEXP threadsSEXP, SEXP accuracySEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< NumericMatrix >::type points(pointsSEXP);
+    Rcpp::traits::input_parameter< NumericMatrix >::type vertices(verticesSEXP);
+    Rcpp::traits::input_parameter< IntegerMatrix >::type faces(facesSEXP);
+    Rcpp::traits::input_parameter< int >::type threads(threadsSEXP);
+    Rcpp::traits::input_parameter< double >::type accuracy(accuracySEXP);
+    rcpp_result_gen = Rcpp::wrap(c_fast_mesh_winding_number(points, vertices, faces, threads, accuracy));
+    return rcpp_result_gen;
+END_RCPP
+}
+// c_mesh_winding_number
+NumericVector c_mesh_winding_number(NumericMatrix points, NumericMatrix vertices, IntegerMatrix faces, int threads);
+RcppExport SEXP _natcpp_c_mesh_winding_number(SEXP pointsSEXP, SEXP verticesSEXP, SEXP facesSEXP, SEXP threadsSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< NumericMatrix >::type points(pointsSEXP);
+    Rcpp::traits::input_parameter< NumericMatrix >::type vertices(verticesSEXP);
+    Rcpp::traits::input_parameter< IntegerMatrix >::type faces(facesSEXP);
+    Rcpp::traits::input_parameter< int >::type threads(threadsSEXP);
+    rcpp_result_gen = Rcpp::wrap(c_mesh_winding_number(points, vertices, faces, threads));
+    return rcpp_result_gen;
+END_RCPP
+}
 // c_ListofMatrixRows
 List c_ListofMatrixRows(const SEXP& object);
 RcppExport SEXP _natcpp_c_ListofMatrixRows(SEXP objectSEXP) {
@@ -173,6 +203,8 @@ static const R_CallMethodDef CallEntries[] = {
     {"_natcpp_c_ijkpos", (DL_FUNC) &_natcpp_c_ijkpos, 5},
     {"_natcpp_c_sub2ind", (DL_FUNC) &_natcpp_c_sub2ind, 2},
     {"_natcpp_c_coords21dindex", (DL_FUNC) &_natcpp_c_coords21dindex, 5},
+    {"_natcpp_c_fast_mesh_winding_number", (DL_FUNC) &_natcpp_c_fast_mesh_winding_number, 5},
+    {"_natcpp_c_mesh_winding_number", (DL_FUNC) &_natcpp_c_mesh_winding_number, 4},
     {"_natcpp_c_ListofMatrixRows", (DL_FUNC) &_natcpp_c_ListofMatrixRows, 1},
     {"_natcpp_c_listlengths", (DL_FUNC) &_natcpp_c_listlengths, 1},
     {"_natcpp_c_topntail", (DL_FUNC) &_natcpp_c_topntail, 1},

src/fast_inside_mesh.cpp

@@ -0,0 +1,79 @@
+// Fast generalised winding number via libigl's "Fast Winding Numbers for Soups
+// and Clouds" (Barill et al. 2018). A bounding-volume hierarchy is built once
+// over the mesh, then each query point is evaluated in O(log F). This is the
+// accelerated back end intended for large point sets on large meshes; the
+// brute-force c_mesh_winding_number() is the O(P*F) reference.
+//
+// We call libigl's per-point primitive (UT_SolidAngle::computeSolidAngle, the
+// same call libigl's batch overload makes internally) directly inside an
+// RcppThread::parallelFor. This keeps the core count controllable via the
+// `threads` argument -- so CRAN's check-farm limit can be respected -- while
+// distributing points across threads ourselves rather than via libigl's own
+// std::thread pool (whose size is fixed by a process-global singleton and so
+// cannot honour a per-call thread count). Routing each query through libigl's
+// templated fast_winding_number(bvh, acc, p) wrapper instead is dramatically
+// slower, so we avoid it.
+//
+// [[Rcpp::depends(RcppEigen)]]
+#include <RcppEigen.h>
+#include <thread>
+#include <RcppThread.h>
+#include "vendor/igl/fast_winding_number.h"
+
+using namespace Rcpp;
+typedef igl::FastWindingNumber::HDK_Sample::UT_Vector3T<float> UTVec3f;
+
+// [[Rcpp::export]]
+NumericVector c_fast_mesh_winding_number(NumericMatrix points,
+                                         NumericMatrix vertices,
+                                         IntegerMatrix faces,
+                                         int threads = 4,
+                                         double accuracy = 2.0) {
+  const int np = points.nrow();
+  const int nv = vertices.nrow();
+  const int nf = faces.nrow();
+  if (points.ncol() != 3)   stop("points must be an Nx3 matrix");
+  if (vertices.ncol() != 3) stop("vertices must be an Nx3 matrix");
+  if (faces.ncol() != 3)    stop("faces must be an Nx3 matrix");
+
+  Eigen::MatrixXd V(nv, 3);
+  for (int i = 0; i < nv; ++i) {
+    V(i, 0) = vertices(i, 0); V(i, 1) = vertices(i, 1); V(i, 2) = vertices(i, 2);
+  }
+  Eigen::MatrixXi F(nf, 3);
+  for (int i = 0; i < nf; ++i) {
+    const int a = faces(i, 0) - 1, b = faces(i, 1) - 1, c = faces(i, 2) - 1;
+    if (a < 0 || b < 0 || c < 0 || a >= nv || b >= nv || c >= nv)
+      stop("faces contains a vertex index outside [1, nrow(vertices)]");
+    F(i, 0) = a; F(i, 1) = b; F(i, 2) = c;
+  }
+
+  // Precompute the BVH once (Taylor expansion order 2, as in libigl examples).
+  igl::FastWindingNumberBVH bvh;
+  igl::fast_winding_number(V, F, 2, bvh);
+
+  // Copy query points into a plain buffer so the worker threads touch no R data
+  // structures and no Eigen expression templates.
+  std::vector<float> Q(static_cast<size_t>(np) * 3);
+  for (int i = 0; i < np; ++i) {
+    Q[3 * i + 0] = static_cast<float>(points(i, 0));
+    Q[3 * i + 1] = static_cast<float>(points(i, 1));
+    Q[3 * i + 2] = static_cast<float>(points(i, 2));
+  }
+
+  NumericVector out(np);
+  double* out_ptr = out.begin();           // NumericVector is not thread-safe
+  const float acc = static_cast<float>(accuracy);
+  const size_t nThreads = (threads > 0) ? static_cast<size_t>(threads)
+                                        : std::thread::hardware_concurrency();
+
+  // computeSolidAngle is a const, read-only query on the shared BVH, so
+  // concurrent single-point evaluation is thread-safe.
+  RcppThread::parallelFor(0, np, [&](int i) {
+    UTVec3f p;
+    p[0] = Q[3 * i + 0]; p[1] = Q[3 * i + 1]; p[2] = Q[3 * i + 2];
+    out_ptr[i] = bvh.ut_solid_angle.computeSolidAngle(p, acc) / (4.0 * igl::PI);
+  }, nThreads);
+
+  return out;
+}