Added vcg_raycaster

CRAN-SUBMISSION

@@ -1,3 +1,3 @@
-Version: 0.1.8
-Date: 2024-09-03 20:22:09 UTC
-SHA: eb45eac5d5d9f6d6ae0c7fff34e2cc7b890eb7f5
+Version: 0.1.9
+Date: 2024-11-08 17:19:07 UTC
+SHA: 8e9e71be96af24ec2444a59a4113abdf02983ec1

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: ravetools
 Type: Package
 Title: Signal and Image Processing Toolbox for Analyzing Intracranial Electroencephalography Data
-Version: 0.1.9
+Version: 0.1.9.9000
 Language: en-US
 Authors@R: c(
     person("Zhengjia", "Wang", email = "[email protected]", role = c("aut", "cre")),

NAMESPACE

@@ -114,6 +114,7 @@ export(shift_array)
 export(surface_path)
 export(vcg_isosurface)
 export(vcg_mesh_volume)
+export(vcg_raycaster)
 export(vcg_smooth_explicit)
 export(vcg_smooth_implicit)
 export(vcg_sphere)

NEWS.md

@@ -1,3 +1,7 @@
+# ravetools 0.2.0
+
+* Added `vcg_raycaster` to find intersection of rays and given mesh object.
+
 # ravetools 0.1.9
 
 * Hot fix by adding `vctrs` to `Suggests` to fix the issue that fails the unit test due to an update in `testthat` package. 

R/RcppExports.R

@@ -613,6 +613,10 @@ vcgDijkstra <- function(vb_, it_, source, maxdist_) {
     .Call(`_ravetools_vcgDijkstra`, vb_, it_, source, maxdist_)
 }
 
+vcgRaycaster <- function(vb_, it_, rayOrigin, rayDirection, maxDistance, bothSides, threads = 1L) {
+    .Call(`_ravetools_vcgRaycaster`, vb_, it_, rayOrigin, rayDirection, maxDistance, bothSides, threads)
+}
+
 # Register entry points for exported C++ functions
 methods::setLoadAction(function(ns) {
     .Call(`_ravetools_RcppExport_registerCCallable`)

R/vcg.R

@@ -482,3 +482,87 @@ vcg_uniform_remesh <- function(
   )
   return(meshintegrity(out))
 }
+
+
+
+#' @title Cast rays to intersect with mesh
+#' @param x surface mesh
+#' @param ray_origin a matrix with 3 rows or a vector of length 3, the positions
+#' of ray origin
+#' @param ray_direction a matrix with 3 rows or a vector of length 3, the
+#' direction of the ray, will be normalized to length 1
+#' @param max_distance positive maximum distance to cast the normalized ray;
+#' default is infinity. Any invalid distances (negative, zero, or \code{NA})
+#' will be interpreted as unset.
+#' @param both_sides whether to inverse the ray (search both positive and
+#' negative ray directions); default is false
+#' @returns A list of ray casting results: whether any intersection is found,
+#' position and face normal of the intersection, distance of the ray, and the
+#' index of the intersecting face (counted from 1)
+#' @examples
+#'
+#' library(ravetools)
+#' sphere <- vcg_sphere(normals = FALSE)
+#' sphere$vb[1:3, ] <- sphere$vb[1:3, ] + c(10, 10, 10)
+#' vcg_raycaster(
+#'   x = sphere,
+#'   ray_origin = array(c(0, 0, 0, 1, 0, 0), c(3, 2)),
+#'   ray_direction = array(c(1, 1, 1, 1, 1, 1), c(3, 2))
+#' )
+#'
+#' @export
+vcg_raycaster <- function(
+    x, ray_origin, ray_direction, max_distance = Inf, both_sides = FALSE) {
+
+  x <- meshintegrity(mesh = x, facecheck = TRUE)
+
+  if(is.matrix(ray_origin)) {
+    ray_origin <- ray_origin[seq_len(3), , drop = FALSE]
+    ray_direction <- ray_direction[seq_len(3), , drop = FALSE]
+  } else {
+    # assuming ray_origin is a vector of 3
+    ray_origin <- matrix(ray_origin[c(1,2,3)], ncol = 1L)
+    ray_direction <- matrix(ray_direction[c(1,2,3)], ncol = 1L)
+  }
+
+  n_rays <- ncol(ray_origin)
+  if(ncol(ray_direction) != n_rays) {
+    stop("`vcg_raycaster`: number of rays is ", n_rays, " according to `ray_origin`. However `ray_direction` is different number of points. Please make sure these two variables have the same number of elements.")
+  }
+
+  # normalize the ray_direction
+  ray_length <- sqrt(colSums(ray_direction^2))
+  zero_length <- is.na(ray_length) | ray_length == 0
+  if(any(zero_length)) {
+    ray_direction[, zero_length] <- 0
+  }
+  ray_direction[, !zero_length] <- ray_direction[, !zero_length] / ray_length[!zero_length]
+
+
+  stopifnot(length(max_distance) == 1)
+  if(is.na(max_distance) || max_distance <= 0) {
+    max_distance <- Inf
+  }
+
+
+  results <- vcgRaycaster(
+    vb_ = x$vb,
+    it_ = x$it - 1L,
+    rayOrigin = ray_origin,
+    rayDirection = ray_direction,
+    maxDistance = max_distance,
+    bothSides = both_sides,
+    threads =
+  )
+
+  list(
+    has_intersection = as.logical(results$hitFlag),
+    intersection = results$intersectPoints,
+    normals = results$intersectNormals,
+    face_index = results$intersectIndex + 1L,
+    distance = results$castDistance,
+    ray_origin = ray_origin,
+    ray_direction = ray_direction
+  )
+
+}

src/RcppExports.cpp

@@ -3607,6 +3607,23 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// vcgRaycaster
+RcppExport SEXP vcgRaycaster(SEXP vb_, SEXP it_, const Rcpp::NumericVector& rayOrigin, const Rcpp::NumericVector& rayDirection, const float& maxDistance, const bool& bothSides, const int& threads);
+RcppExport SEXP _ravetools_vcgRaycaster(SEXP vb_SEXP, SEXP it_SEXP, SEXP rayOriginSEXP, SEXP rayDirectionSEXP, SEXP maxDistanceSEXP, SEXP bothSidesSEXP, SEXP threadsSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< SEXP >::type vb_(vb_SEXP);
+    Rcpp::traits::input_parameter< SEXP >::type it_(it_SEXP);
+    Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type rayOrigin(rayOriginSEXP);
+    Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type rayDirection(rayDirectionSEXP);
+    Rcpp::traits::input_parameter< const float& >::type maxDistance(maxDistanceSEXP);
+    Rcpp::traits::input_parameter< const bool& >::type bothSides(bothSidesSEXP);
+    Rcpp::traits::input_parameter< const int& >::type threads(threadsSEXP);
+    rcpp_result_gen = Rcpp::wrap(vcgRaycaster(vb_, it_, rayOrigin, rayDirection, maxDistance, bothSides, threads));
+    return rcpp_result_gen;
+END_RCPP
+}
 
 // validate (ensure exported C++ functions exist before calling them)
 static int _ravetools_RcppExport_validate(const char* sig) { 
@@ -3930,6 +3947,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_ravetools_vcgVolume", (DL_FUNC) &_ravetools_vcgVolume, 1},
     {"_ravetools_vcgSphere", (DL_FUNC) &_ravetools_vcgSphere, 2},
     {"_ravetools_vcgDijkstra", (DL_FUNC) &_ravetools_vcgDijkstra, 4},
+    {"_ravetools_vcgRaycaster", (DL_FUNC) &_ravetools_vcgRaycaster, 7},
     {"_ravetools_RcppExport_registerCCallable", (DL_FUNC) &_ravetools_RcppExport_registerCCallable, 0},
     {NULL, NULL, 0}
 };

src/vcgCommon.cpp

@@ -568,3 +568,164 @@ SEXP vcgDijkstra(SEXP vb_, SEXP it_, const Rcpp::IntegerVector & source, const d
   }
   return R_NilValue; // -Wall
 }
+
+
+// [[Rcpp::export]]
+RcppExport SEXP vcgRaycaster(
+    SEXP vb_ , SEXP it_,
+    const Rcpp::NumericVector & rayOrigin, // 3 x n matrix
+    const Rcpp::NumericVector & rayDirection,
+    const float & maxDistance,
+    const bool & bothSides,
+    const int & threads = 1)
+{
+  try {
+    ravetools::MyMesh m;
+    int check = ravetools::IOMesh<ravetools::MyMesh>::vcgReadR(m, vb_, it_);
+    if (check != 0) {
+      throw std::runtime_error("Mesh has no faces or vertices. Unable to perform raycaster");
+    }
+
+    ravetools::ScalarType x,y,z;
+    ravetools::MyMesh rays;
+
+    // Leave R to check if nRays > 0
+    unsigned int nRays = rayOrigin.length() / 3;
+    Rcpp::NumericVector castDistance(nRays);
+    Rcpp::IntegerVector hitFlag(nRays);
+
+    // rayOrigin and rayDirection must be identical 3xn
+    // Leave the checks in R wrapper
+    Rcpp::NumericVector intersectPoints(rayOrigin);
+    Rcpp::NumericVector intersectNormals(rayDirection);
+    Rcpp::IntegerVector intersectIndex(nRays);
+
+    //Allocate target
+    std::vector<ravetools::MyMesh::VertexPointer> ivp;
+    vcg::tri::Allocator<ravetools::MyMesh>::AddVertices(rays, nRays);
+    vcg::Point3f normtmp;
+
+    // Copy the rayOrigin and rayDirection
+    Rcpp::NumericVector::iterator intersectPointsPtr = intersectPoints.begin();
+    Rcpp::NumericVector::iterator intersectNormalsPtr = intersectNormals.begin();
+    ravetools::MyMesh::VertexIterator vi = rays.vert.begin();
+    for (unsigned int i=0; i < nRays; i++, vi++) {
+      x = *intersectPointsPtr++;
+      y = *intersectPointsPtr++;
+      z = *intersectPointsPtr++;
+      (*vi).P() = ravetools::MyMesh::CoordType(x, y, z);
+      x = *intersectNormalsPtr++;
+      y = *intersectNormalsPtr++;
+      z = *intersectNormalsPtr++;
+      (*vi).N() = ravetools::MyMesh::CoordType(x, y, z);
+    }
+
+    // bounding box to calculate max cast distance
+    vcg::tri::UpdateBounding<ravetools::MyMesh>::Box(m);
+    m.face.EnableNormal();
+    vcg::tri::UpdateNormal<ravetools::MyMesh>::PerFaceNormalized(m);
+    vcg::tri::UpdateNormal<ravetools::MyMesh>::PerVertexAngleWeighted(m);
+    vcg::tri::UpdateNormal<ravetools::MyMesh>::NormalizePerVertex(m);
+
+    vcg::tri::UpdateNormal<ravetools::MyMesh>::NormalizePerVertex(rays);
+
+    vcg::tri::FaceTmark<ravetools::MyMesh> mf;
+    mf.SetMesh( &m );
+    vcg::RayTriangleIntersectionFunctor<true> FintFunct;
+    vcg::GridStaticPtr<ravetools::MyMesh::FaceType, ravetools::MyMesh::ScalarType> gridSearcher;
+    gridSearcher.Set(m.face.begin(), m.face.end());
+
+#pragma omp parallel for firstprivate(maxDistance,gridSearcher,mf) schedule(static) num_threads(threads)
+{
+    for ( int i = 0; i < rays.vn ; i++ ) {
+      float t0 = 0.0f, t1 = 0.0f;
+      int faceIndex = -1;
+      vcg::Ray3f ray;
+      vcg::Point3f orig = rays.vert[i].P();
+      // vcg::Point3f orig0 = orig;
+      vcg::Point3f dir = rays.vert[i].N();
+      vcg::Point3f intersection = ravetools::MyMesh::CoordType(0, 0, 0);
+      ravetools::MyFace* facePtr0;
+      ravetools::MyFace* facePtr1;
+
+      /**
+       *  Set ray origin to be slightly "behind" the `orig`
+       *  This is because if orig coincide with the underlying intersection,
+       *  FintFunct will not be able to identify the intersection
+       *  Example:
+       sphere <- ravetools::vcg_sphere()
+       box <- Rvcg::vcgBox(sphere)
+       box$vb[1:3,] <- box$vb[1:3,] + c(1,1,1) - 1e-6
+       mesh <- box
+       vcgRaycaster(vb_ = mesh$vb, it_ = mesh$it - 1L, rayOrigin = matrix(c(0,0,0), ncol = 1), rayDirection = matrix(c(1,1,1), ncol = 1), maxDistance = 1e14, bothSides = FALSE)
+       */
+      ray.SetOrigin(orig - 1e-6f * dir);
+      ray.SetDirection(dir);
+
+      // raycaster
+      facePtr0 = GridDoRay(gridSearcher, FintFunct, mf, ray, maxDistance, t0);
+      if ( bothSides ) {
+        ray.SetOrigin(orig + 1e-6f * dir);
+        // cast the ray backwards
+        ray.SetDirection(-dir);
+        facePtr1 = GridDoRay(gridSearcher, FintFunct, mf, ray, maxDistance, t1);
+        if( facePtr1 && ( !facePtr0 || t1 < t0 ) ) {
+          facePtr0 = facePtr1;
+          t0 = -t1;
+        }
+      }
+
+      if( facePtr0 ) {
+        // pay off the debt
+        if( t0 > 0.0f ) {
+          t0 -= 1e-6f;
+        } else {
+          t0 += 1e-6f;
+        }
+
+        intersection = rays.vert[i].P()+dir * t0;
+        castDistance[ i ] = t0;
+        hitFlag[ i ] = 1;
+
+        faceIndex = vcg::tri::Index(m, facePtr0);
+
+        // face normal
+        const ravetools::MyFace face = m.face[faceIndex];
+        ravetools::MyMesh::CoordType faceNormal = (face.V(0)->N() + face.V(1)->N() + face.V(2)->N()).normalized();
+
+        intersectPoints[i * 3] = intersection[0];
+        intersectPoints[i * 3 + 1] = intersection[1];
+        intersectPoints[i * 3 + 2] = intersection[2];
+        intersectNormals[i * 3] = faceNormal[0];
+        intersectNormals[i * 3 + 1] = faceNormal[1];
+        intersectNormals[i * 3 + 2] = faceNormal[2];
+        intersectIndex[i] = faceIndex;
+
+      } else {
+        // No intersection
+        castDistance[ i ] = NA_REAL;
+        hitFlag[ i ] = 0;
+
+        intersectPoints[i * 3] = NA_REAL;
+        intersectPoints[i * 3 + 1] = NA_REAL;
+        intersectPoints[i * 3 + 2] = NA_REAL;
+        intersectNormals[i * 3] = NA_REAL;
+        intersectNormals[i * 3 + 1] = NA_REAL;
+        intersectNormals[i * 3 + 2] = NA_REAL;
+        intersectIndex[i] = NA_INTEGER;
+      }
+    }
+}
+    return Rcpp::List::create(Rcpp::Named("intersectPoints") = intersectPoints,
+                              Rcpp::Named("intersectNormals") = intersectNormals,
+                              Rcpp::Named("intersectIndex") = intersectIndex,
+                              Rcpp::Named("hitFlag") = hitFlag,
+                              Rcpp::Named("castDistance") = castDistance
+    );
+  } catch (std::exception& e) {
+    Rcpp::stop( e.what());
+  } catch (...) {
+    Rcpp::stop("unknown exception");
+  }
+  return R_NilValue;
+}