In Biology and Physical Anthropology, sliding semi-landmarks are widely used to estimate the surface shape of a biological structure. The...

KD-trees are a fast and reliable method for spatial indexing and closest point searches. In Rvcg, the functions vcgKDtree and vcgClostKD provide these functionality, however only returning the result and discarding the KD-trees. For registration processes, the distances are updated in each iteration and it would be sensible to save and reuse the computationally expensive KD-tree setup, especially for large meshes.

How to achieve this

The package Rcpp, that is used to exchange data between R and vcglib, is also sporting external pointers (some sort of smart pointers) that allow to return pointers to objects created by the C++ code and later pass it back and reuse it. The newly introduced function vcgCreateKD tree allows to create such an object from meshes and point clouds. Using vcgSearchKDtree allows to calculated the k-closest neighbours and the respective distances. For closest point searches on triangular meshes, this is based on visiting the faces with the k-closest barycenters. Apart from the KD-tree and the corresponding points (the faces’ barycenters in this case), we now also need the information of the surface mesh. This can be obtained by the function vcgCreateKDtreeFromBarycenters and to search it using the function vcgClostOnKDtreeFromBarycenters can be used.

Below is the speed comparison of 10 closest point searches for 5637 coordinates on a triangular mesh with 564877 vertices and 1127080 faces.

CAVEAT: In order to reproduce this, you need to install the xptr branch from Rvcg: devtools::install_github("zarquon42b/Rvcg",ref="xptr")

Here the results (Figure 1):

• setting up the kd-tree takes 1.7 secs
• running 10 searches without reusing the KD-tree takes 11.3 secs
• running 10 searches reusing the KD-tree takes 11.3 secs only 3.7 secs

require(Rvcg)
data(humface)
data(dummyhead)
## do a face subdivision to get a mesh with 564877 vertices and 1127080 faces
humfacehigh <- vcgSubdivide(humface,threshold = 0.5,iterations = 8)

tbary <- system.time(kdtreeBary <- vcgCreateKDtreeFromBarycenters(humfacehigh))
print(tbary[3])


tsearch0 <- system.time(lapply(1:10,function(x) test <- vcgClostOnKDtreeFromBarycenters(kdtreeBary,dummyhead.mesh,threads=parallel::detectCores())))
print(tsearch0[3])

## now without reusing the kd-tree
tsearch1 <- system.time(lapply(1:10,function(x) test <- vcgClostKD(dummyhead.mesh,humfacehigh,threads=parallel::detectCores())))
print(tsearch1[3])

General stuff

As can be seen in my last four posts, I started playing with displacement fields and interpolation....

As RvtkStatismo is already sporting a B-spline kernel for Gaussian-Process models (a modified version contained in the statismo command...

To make the plotting of displacement fields more comprehensive, I deprecated plotDisplacementField in favor of a plotDisplacementField in favor of a plot-method for objects of class DisplacementField, returning (and optionally plotting) an object of the (newly introduced) class DisplacementPlot, that in turn has its own plot-method.

For better visualization, and as rgl has no rendering options for arrows, I decided to render the starting positions as points, leading to lollipop graphs (Fig. 1 shows the resulting plot, using the example from yesterday). The code for the plot below is simply:

plot(dispfield)
wire3d(dummytrans)

And finally, there is a new complementing function invertDisplacementField to invert an existing displacement field.