obbox & sampleInMesh

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: cgalMeshes
 Title: R6 Based Utilities for 3D Meshes using 'CGAL'
-Version: 2.2.0.9001
+Version: 2.2.0.9002
 Authors@R: 
  person("Stéphane", "Laurent", , "[email protected]", role = c("aut", "cre"))
 Maintainer: Stéphane Laurent <[email protected]>

NEWS.md

@@ -4,6 +4,10 @@
 
 - Implementation of the 3D alpha wrapping.
 
+- Optimal bounding box of a mesh.
+
+- New method `$sampleInMesh` to uniformly sample in the volume bounded by a mesh. The `$sample` method has been renamed to `$sampleOnMesh`.
+
 
 # cgalMeshes 2.2.0
 

R/cgalMesh.R

@@ -1339,7 +1339,26 @@ cgalMesh <- R6Class(
  . <- private[[".CGALmesh"]]$merge(mesh2XPtr)
  invisible(self)
  },
-
+
+ #' @description Optimal bounding box.
+ #' @return A \strong{rgl} mesh of the optimal hexahedron enclosing the mesh. 
+ #' @examples 
+ #' \donttest{library(rgl)
+ #' cyclide <- cyclideMesh(a = 97, c = 32, mu = 57)
+ #' mesh <- cgalMesh$new(cyclide)
+ #' obb <- mesh$optimalBoundingBox()
+ #' open3d(windowRect = 50 + c(0, 0, 512, 512), zoom = 0.9)
+ #' shade3d(cyclide, color = "green")
+ #' shade3d(obb, color = "yellow", alpha = 0.25)}
+ "optimalBoundingBox" = function() {
+ hx <- private[[".CGALmesh"]]$optimalBoundingBox()
+ obb <- hx[["rmesh"]]
+ qmesh3d(
+ vertices = obb[["vertices"]],
+ indices = obb[["faces"]]
+ )
+ },
+
  #' @description Reorient the connected components of the mesh in order that 
  #' it bounds a volume. The mesh must be triangle.
  #' @return The modified \code{cgalMesh} object, invisibly. \strong{WARNING}: 
@@ -1399,14 +1418,27 @@ cgalMesh <- R6Class(
  invisible(self)
  },
 
+ #' @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; 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.
- "sample" = function(nsims) {
+ "sampleOnMesh" = function(nsims) {
  stopifnot(isStrictPositiveInteger(nsims))
- private[[".CGALmesh"]]$sampleMesh(as.integer(nsims))
+ private[[".CGALmesh"]]$sampleOnMesh(as.integer(nsims))
  },
-
+  
  #' @description Check whether the mesh self-intersects. The mesh must be 
  #' triangle.
  #' @return A Boolean value, whether the mesh self-intersects.

src/MODULE.cpp

@@ -68,11 +68,13 @@ RCPP_MODULE(class_CGALmesh) {
  .method("isValidPolygonMesh", &CGALmesh::isValidPolygonMesh)
  .method("LoopSubdivision", &CGALmesh::LoopSubdivision)
  .method("merge", &CGALmesh::merge)
+ .method("optimalBoundingBox", &CGALmesh::optimalBoundingBox)
  .method("orientToBoundVolume", &CGALmesh::orientToBoundVolume)
  .method("print", &CGALmesh::print)
  .method("removeSelfIntersections", &CGALmesh::removeSelfIntersections)
  .method("reverseFaceOrientations", &CGALmesh::reverseFaceOrientations)
- .method("sampleMesh", &CGALmesh::sampleMesh)
+ .method("sampleInMesh", &CGALmesh::sampleInMesh)
+ .method("sampleOnMesh", &CGALmesh::sampleOnMesh)
  .method("sharpEdges", &CGALmesh::sharpEdges)
  .method("Sqrt3Subdivision", &CGALmesh::Sqrt3Subdivision)
  .method("subtract", &CGALmesh::subtract)

src/MODULE.h

@@ -1648,6 +1648,37 @@ class CGALmesh {
  }
 
 
+ // ----------------------------------------------------------------------- //
+ // ----------------------------------------------------------------------- //
+ Rcpp::List optimalBoundingBox() {
+ Mesh3 kmesh;
+ CGAL::copy_face_graph(mesh, kmesh);
+ std::array<Point3, 8> obb_points;
+ CGAL::oriented_bounding_box(kmesh, obb_points,
+ CGAL::parameters::use_convex_hull(true));
+ // Make a mesh out of the oriented bounding box
+ Mesh3 obbMesh;
+ CGAL::make_hexahedron(
+ obb_points[0], obb_points[1], obb_points[2], obb_points[3],
+ obb_points[4], obb_points[5], obb_points[6], obb_points[7], obbMesh
+ );
+ EMesh3 obbEMesh;
+ CGAL::copy_face_graph(obbMesh, obbEMesh);
+ Rcpp::List rmesh = RSurfEKMesh2(obbEMesh, false, 4);
+ Rcpp::NumericMatrix hxVertices(3, 8);
+ for(int i = 0; i < 8; i++) {
+ Point3 pt = obb_points[i];
+ Rcpp::NumericVector v = 
+ Rcpp::NumericVector::create(pt.x(), pt.y(), pt.z());
+ hxVertices(Rcpp::_, i) = v;
+ }
+ return Rcpp::List::create(
+ Rcpp::Named("rmesh") = rmesh,
+ Rcpp::Named("hxVertices") = hxVertices
+ );
+ }
+
+
  // ----------------------------------------------------------------------- //
  void orientToBoundVolume() {
  if(!CGAL::is_triangle_mesh(mesh)) {
@@ -1689,7 +1720,73 @@ class CGALmesh {
 
  // ----------------------------------------------------------------------- //
  // ----------------------------------------------------------------------- //
- Rcpp::NumericMatrix sampleMesh(const unsigned nsims) {
+ Rcpp::NumericMatrix sampleInMesh(const unsigned nsims) {
+
+ if(!CGAL::is_triangle_mesh(mesh)) {
+ Rcpp::stop("The mesh is not triangle.");
+ }
+ if(!CGAL::is_closed(mesh)) {
+ Rcpp::stop("The mesh is not closed.");
+ }
+
+ Mesh3 kmesh = epeck2epick(mesh);
+
+ CGAL::Side_of_triangle_mesh<Mesh3, K> Where(kmesh);
+
+ std::array<Point3, 8> hxh;
+ CGAL::oriented_bounding_box(
+ kmesh, hxh, CGAL::parameters::use_convex_hull(true)
+ );
+
+ std::array<std::array<Vector3, 4>, 5> ths = hexahedronTetrahedra(hxh);
+ std::array<Vector3, 4> th1 = ths[0];
+ std::array<Vector3, 4> th2 = ths[1];
+ std::array<Vector3, 4> th3 = ths[2];
+ std::array<Vector3, 4> th4 = ths[3];
+ std::array<Vector3, 4> th5 = ths[4];
+ const double vol1 = volumeTetrahedron(
+ V3toP3(th1[0]), V3toP3(th1[1]), V3toP3(th1[2]), V3toP3(th1[3])
+ );
+ const double vol2 = volumeTetrahedron(
+ V3toP3(th2[0]), V3toP3(th2[1]), V3toP3(th2[2]), V3toP3(th2[3])
+ );
+ const double vol3 = volumeTetrahedron(
+ V3toP3(th3[0]), V3toP3(th3[1]), V3toP3(th3[2]), V3toP3(th3[3])
+ );
+ const double vol4 = volumeTetrahedron(
+ V3toP3(th4[0]), V3toP3(th4[1]), V3toP3(th4[2]), V3toP3(th4[3])
+ );
+ const double vol5 = volumeTetrahedron(
+ V3toP3(th5[0]), V3toP3(th5[1]), V3toP3(th5[2]), V3toP3(th5[3])
+ );
+ Rcpp::NumericVector volumes = Rcpp::NumericVector::create(
+ vol1, vol2, vol3, vol4, vol5
+ );
+ 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) {
+ int index = die5(gen);
+ std::array<Vector3, 4> th = ths[index];
+ Vector3 v = sampleTetrahedron(th[0], th[1], th[2], th[3], gen);
+ Point3 p = V3toP3(v);
+ CGAL::Bounded_side side = Where(p);
+ if(side == CGAL::ON_BOUNDED_SIDE || side == CGAL::ON_BOUNDARY) {
+ Rcpp::NumericVector sim = 
+ Rcpp::NumericVector::create(p.x(), p.y(), p.z());
+ Sims(Rcpp::_, i++) = sim;
+ }
+ }
+ return Rcpp::transpose(Sims); 
+ } 
+
+
+ // ----------------------------------------------------------------------- //
+ // ----------------------------------------------------------------------- //
+ Rcpp::NumericMatrix sampleOnMesh(const unsigned nsims) {
  if(!CGAL::is_triangle_mesh(mesh)) {
  Rcpp::stop("The mesh is not triangle.");
  }

src/cgalMesh.h

@@ -79,6 +79,12 @@
 #include <CGAL/Side_of_triangle_mesh.h>
 #include <CGAL/alpha_wrap_3.h>
 
+#include <CGAL/optimal_bounding_box.h>
+
+#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/discrete_distribution.hpp>
+#include <boost/random/uniform_real_distribution.hpp>
+
 // -------------------------------------------------------------------------- //
 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
 typedef K::Point_3 Point3;
@@ -215,6 +221,14 @@ void copy_property(
 Mesh3 epeck2epick(EMesh3&);
 EMesh3 epick2epeck(Mesh3&);
 
+double volumeTetrahedron(Point3, Point3, Point3, Point3);
+Vector3 P3toV3(Point3);
+Point3 V3toP3(Vector3);
+std::array<std::array<Vector3, 4>, 5> hexahedronTetrahedra(
+ std::array<Point3, 8>
+);
+Vector3 sampleTetrahedron(Vector3, Vector3, Vector3, Vector3, boost::mt19937);
+
 
 ////////////////////////////////////////////////////////////////////////////////