Splitting RCGAL… and the connected components
Want to share your content on Rbloggers? click here if you have a blog, or here if you don't.
I published two posts here about my package RCGAL, the first one about surface reconstruction, and the second one about Boolean operations on 3D meshes.
Now I decided to split this package into two packages: SurfaceReconstruction, and MeshesOperations.
I think MeshesOperations is almost ready for submission to CRAN. I have just started SurfaceReconstruction today, but it already works.
Dealing with a ‘R CMD check’ issue
A technical note before playing with the MeshesOperations package. Feel free to skip this section.
I faced a ‘R CMD check’ issue on ‘winbuilder’ with MeshesOperations. Namely, this socalled significant warning:
C:/.../BH/include/boost/container/detail/copy_move_algo.hpp:184:19: warning: 'void* memmove(void*, const void*, size_t)' writing to an object of type 'value_type' {aka 'struct std::pair, std::allocator >, long long unsigned int> > >, std::allocator, std::allocator >, long long unsigned int> > > > >, CGAL::internal::In_place_list_iterator, std::allocator >, long long unsigned int> > >, std::allocator, std::allocator >, long long unsigned int> > > > > >'} with no trivial copyassignment; use copyassignment or copyinitialization instead [Wclassmemaccess]
This is a warning from the BH package (which allows to use the C++ library Boost with Rcpp). On Linux, the ‘R CMD check’ report was clean.
I found a solution to this issue. I describe it here, in case it could help someone.
In order to use the BH package, I included it in the LinkingTo field of the DESCRIPTION file. Then I discovered that this was not necessary with R4.2.0, and that’s because Rtools42 contains the Boost library.
So this solved the problem for Windows, but created a new problem for the other OSs. Then I solved this new problem by:

putting BH in the Suggests field of DESCRIPTION (and only in this field)

adding the following lines in the Makevars file in the src folder of the package (breaking the first line for display here):
BH_PATH = `echo 'cat(system.file("include", package = "BH", mustWork=TRUE))' \  "${R_HOME}/bin/R" vanilla noecho` PKG_CXXFLAGS = I$(BH_PATH)
The first line allows to call R to get the path of the include folder of the BH package, and this path is set to the variable BH_PATH.
Of course I didn’t change the Makevars.win file. This works, but I don’t know yet whether this is acceptable for CRAN. However I am optimistic, since such a technique can be found in ‘Writing R extensions’.
Connected components
There’s something new in MeshesOperations as compared to RCGAL: the computation of the connected components of a mesh. I’m going to show you an example. Note that the mesh in this example is triangle (i.e. it only has triangular faces), and then one could alternatively use my package MeshesOperations (not submitted to CRAN yet) to get the connected components of this mesh.
So let’s try this feature. Our example is an isosurface called the C8 surface. Unfortunately I do not remember where I found it. I’m using the rmarchingcubes to construct a mesh of this isosurface.
f < function(x, y, z){ 64*x**8  128*x**6 + 80*x**4  16*x**2 + 2 + 64*y**8  128*y**6 + 80*y**4  16*y**2 + 64*z**8  128*z**6 + 80*z**4  16*z**2 } ngrid < 150L x < y < z < seq(1.1, 1.1, len = ngrid) Grid < expand.grid(X = x, Y = y, Z = z) voxel < array( with(Grid, f(X, Y, Z)), dim = c(ngrid, ngrid, ngrid) ) library(rmarchingcubes) contour_shape < contour3d( griddata = voxel, level = 0.1, x = x, y = y, z = z )
Now let’s plot it with rgl:
library(rgl) tmesh < tmesh3d( vertices = t(contour_shape[["vertices"]]), indices = t(contour_shape[["triangles"]]), normals = contour_shape[["normals"]], homogeneous = FALSE ) open3d(windowRect = c(50, 50, 562, 562)) shade3d(tmesh, color = "darkred")
As you can see, it has many isolated parts, the socalled connected components.
Now we use the MeshesOperations package to extract these connected components, while requesting the vertex normals of each component:
library(MeshesOperations) meshes < connectedComponents( vertices = contour_shape[["vertices"]], faces = contour_shape[["triangles"]], normals = TRUE ) # Found 64 components.
Now let’s plot all theses meshes, one color per mesh:
ncc < length(meshes) library(randomcoloR) colors < randomColor(ncc, hue = "random", luminosity = "dark") open3d(windowRect = c(50, 50, 562, 562)) for(i in 1L:ncc){ cc < meshes[[i]] tmesh < tmesh3d( vertices = t(cc[["vertices"]]), indices = t(cc[["faces"]]), normals = cc[["normals"]], homogeneous = FALSE ) shade3d(tmesh, color = colors[i]) }
Nice. However there is one problem: compare the smoothness of these meshes to the smoothness of our first plot. It is not as smooth.
I think I know why. The smoothness depends on the vertex normals. Here
we computed the normals with CGAL, then each vertex normal is obtained
by averaging the normals of the surrounding faces of this vertex. We
could have used rgl::addNormals
instead, the result would
be the same. But there is a convenient way to get the “true” vertex
normals of an isosurface: the normal associated to a vertex is the
gradient of the isosurface function evaluated at the coordinates of this
vertex. So my guess is that rmarchingcubes uses the
gradient. I am still thinking of having the possibility to preserve the
original normals when computing the connected components, I didn’t have
the time yet to find how to do that. This is not a problem, we will
compute the true normals with R. In fact we will probably do better than
rmarchingcubes: surely this package numerically
computes the gradient, while we will use the exact gradient.
How? Look at the isosurface function f
: this is a
multivariate polynomial. So it is easy to calculate its derivatives.
This would be tedious to calculate them by hand however, so let’s use
the spray package:
library(spray) # define the polynomial corresponding to f: P < f(lone(1,3), lone(2,3), lone(3,3)) # and its derivatives: dfx < as.function(deriv(P, 1L)) dfy < as.function(deriv(P, 2L)) dfz < as.function(deriv(P, 3L)) # the gradient is: gradient < function(xyz){ cbind(dfx(xyz), dfy(xyz), dfz(xyz)) }
Just a couple of lines of code… Now let’s do the plot with the true normals:
open3d(windowRect = c(50, 50, 562, 562)) for(i in 1:ncc){ cc < meshes[[i]] tmesh < tmesh3d( vertices = t(cc$vertices), indices = t(cc$faces), normals = gradient(cc$vertices), homogeneous = FALSE ) shade3d(tmesh, color = colors[i]) }
A perfect smoothness!
Rbloggers.com offers daily email updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/datascience job.
Want to share your content on Rbloggers? click here if you have a blog, or here if you don't.