---
title: "Splitting RCGAL... and the connected components"
author: "Stéphane Laurent"
date: '2022-05-27'
tags: geometry, R, rgl, graphics, C++, Rcpp
rbloggers: yes
output:
md_document:
variant: markdown
preserve_yaml: true
html_document:
highlight: kate
keep_md: no
highlighter: pandoc-solarized
---
I published two posts here about my package **RCGAL**, the
[first one](https://laustep.github.io/stlahblog/posts/SurfaceReconstruction.html)
about surface reconstruction, and the
[second one](https://laustep.github.io/stlahblog/posts/BooleanOpsOnMeshes.html)
about Boolean operations on 3D meshes.
Now I decided to split this package into two packages:
[SurfaceReconstruction](https://github.com/stla/SurfaceReconstruction),
and
[MeshesOperations](https://github.com/stla/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 'win-builder' with **MeshesOperations**.
Namely, this so-called 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 copy-assignment; use copy-assignment or copy-initialization
instead [-Wclass-memaccess]
```
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 R-4.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 --no-echo`
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](https://github.com/stla/concomp) (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.
```r
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:
```r
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")
```
![](./figures/C8surface.gif)
As you can see, it has many isolated parts, the so-called
*connected components*.
Now we use the **MeshesOperations** package to extract these connected
components, while requesting the vertex normals of each component:
```r
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:
```r
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])
}
```
![](./figures/C8_components.gif)
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:
```r
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:
```r
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])
}
```
![](./figures/C8_components_trueNormals.gif)
A perfect smoothness!