cleaning

R/cgalMesh.R

@@ -1443,24 +1443,26 @@ cgalMesh <- R6Class(
       invisible(self)
     },
 
-    #' @description Random sampling on the mesh. The mesh must be triangle.
-    #' @param nsims integer, the desired number of simulations
-    #' @return A \code{nsims x 3} matrix containing the simulations.
-    "sampleOnMesh" = function(nsims) {
-      stopifnot(isStrictPositiveInteger(nsims))
-      private[[".CGALmesh"]]$sampleOnMesh(as.integer(nsims))
-    },
-
     #' @description Random sampling in the volume bounded by the mesh. The 
     #'   mesh must be closed and triangle. The method consists in sampling in 
     #'   the optimal bounding box of the mesh and rejecting the points that 
-    #'   do not fall inside the volume.
+    #'   do not fall inside the volume; therefore this is inefficient when 
+    #'   the volume of the mesh is small as compared to the volume of its 
+    #'   bounding box.
     #' @param nsims integer, the desired number of simulations
     #' @return A \code{nsims x 3} matrix containing the simulations.
     "sampleInMesh" = function(nsims) {
       stopifnot(isStrictPositiveInteger(nsims))
       private[[".CGALmesh"]]$sampleInMesh(as.integer(nsims))
     },
+
+    #' @description Random sampling on the mesh. The mesh must be triangle.
+    #' @param nsims integer, the desired number of simulations
+    #' @return A \code{nsims x 3} matrix containing the simulations.
+    "sampleOnMesh" = function(nsims) {
+      stopifnot(isStrictPositiveInteger(nsims))
+      private[[".CGALmesh"]]$sampleOnMesh(as.integer(nsims))
+    },
 
     #' @description Check whether the mesh self-intersects. The mesh must be 
     #'   triangle.

src/MODULE.h

@@ -1986,7 +1986,7 @@ class CGALmesh {
   // ----------------------------------------------------------------------- //
   // ----------------------------------------------------------------------- //
   Rcpp::NumericMatrix sampleInMesh(const unsigned nsims) {
-    boost::mt19937 gen;
+
     if(!CGAL::is_triangle_mesh(mesh)) {
       Rcpp::stop("The mesh is not triangle.");
     }
@@ -2030,6 +2030,7 @@ class CGALmesh {
     Rcpp::NumericVector probs = volumes / sum(volumes);
     boost::random::discrete_distribution<> die5(probs.begin(), probs.end());
     // sampling
+    boost::mt19937 gen;
     Rcpp::NumericMatrix Sims(3, nsims);
     unsigned i = 0;
     while(i < nsims) {

src/unexported_sampling.cpp

@@ -3,23 +3,34 @@
 #endif
 
 
+// ------------------------------------------------------------------------- //
+// - volume of tetrahedron p1p2p3p4 ---------------------------------------- //
 double volumeTetrahedron(Point3 p1, Point3 p2, Point3 p3, Point3 p4) {
   Vector3 v = p1 - p4;
   Vector3 w = CGAL::cross_product(p2 - p4, p3 - p4);
   double dotprod = CGAL::scalar_product(v, w);
   return fabs(dotprod) / 6.0;
 }
 
+
+// ------------------------------------------------------------------------- //
+// - point to vector ------------------------------------------------------- //
 Vector3 P3toV3(Point3 p) {
   Vector3 v(p.x(), p.y(), p.z());
   return v;
 }
 
+
+// ------------------------------------------------------------------------- //
+// - vector to point ------------------------------------------------------- //
 Point3 V3toP3(Vector3 v) {
   Point3 p(v.x(), v.y(), v.z());
   return p;
 }
 
+
+// ------------------------------------------------------------------------- //
+// - five tetrahedra partitioning a hexahedron ----------------------------- //
 std::array<std::array<Vector3, 4>, 5> hexahedronTetrahedra(
     std::array<Point3, 8> hxh
 ) {
@@ -40,6 +51,9 @@ std::array<std::array<Vector3, 4>, 5> hexahedronTetrahedra(
   return tetrahedra;
 }
 
+
+// ------------------------------------------------------------------------- //
+// - sample one point in tetrahedron v1v2v3v4  ----------------------------- //
 Vector3 sampleTetrahedron(
   Vector3 v1, Vector3 v2, Vector3 v3, Vector3 v4, boost::mt19937 gen
 ) {