3d surface plot with xyz coordinates

Here is an approach using kernel density estimation and the contour3d function from misc3d. I played around until I found a value for levels that worked decently. It's not perfectly precise, but you may be able to tweak things to get a better, more accurate surface. If you have more than 8GB of memory, then you may be able to increase n beyond what I did here.

library(rgl)library(misc3d)library(onion); data(bunny)# the larger the n, the longer it takes, the more RAM you needbunny.dens <- kde3d(bunny[,1],bunny[,2],bunny[,3], n=150,     lims=c(-.1,.2,-.1,.2,-.1,.2)) # I chose lim values manuallycontour3d(bunny.dens$d, level = 600,     color = "pink", color2 = "green", smooth=500)rgl.viewpoint(zoom=.75)

The image on the right is from the bottom, just to show another view.

You can use a larger value for n in kde3d but it will take longer, and you may run out of RAM if the array becomes too large. You could also try a different bandwidth (default used here). I took this approach from Computing and Displaying Isosurfaces in R - Feng & Tierney 2008.

Very similar isosurface approach using the Rvcg package:

library(Rvcg)library(rgl)library(misc3d)library(onion); data(bunny)bunny.dens <- kde3d(bunny[,1],bunny[,2],bunny[,3], n=150,     lims=c(-.1,.2,-.1,.2,-.1,.2)) # I chose lim values manuallybunny.mesh <- vcgIsosurface(bunny.dens$d, threshold=600)shade3d(vcgSmooth(bunny.mesh,"HC",iteration=3),col="pink") # do a little smoothing

Since it's a density estimation based approach, we can get a little more out of it by increasing the density of the bunny. I also use n=400 here. The cost is a significant increase in computation time, but the resulting surface is a hare better:

bunny.dens <- kde3d(rep(bunny[,1], 10), # increase density.                    rep(bunny[,2], 10),                    rep(bunny[,3], 10), n=400,                     lims=c(-.1,.2,-.1,.2,-.1,.2))bunny.mesh <- vcgIsosurface(bunny.dens$d, threshold=600)shade3d(vcgSmooth(bunny.mesh,"HC",iteration=1), col="pink")

Better, more efficient surface reconstruction methods exist (e.g. power crust, Poisson surface reconstruction, ball-pivot algorithm), but I don't know that any have been implemented in R, yet.

Here's a relevant Stack Overflow post with some great information and links to check out (including links to code): robust algorithm for surface reconstruction from 3D point cloud?.

I think have found one possible solution using the alphashape3d package. I had to play around a bit to get an acceptable value for alpha, which is related to the distances in the given data set (e.g. sd of bunny gave me some insight). I'm still trying to figure out how to better control the width of lines in vertices and edges so as not to dominate the plot, but this is probably related to settings in rgl.

Example:

library(onion)library(rgl)library(geometry)library(alphashape3d)data(bunny)apply(bunny,2,sd)alphabunny <- ashape3d(bunny, alpha = 0.003)bg3d(1)plot.ashape3d(alphabunny, col=c(5,5,5), lwd=0.001, size=0, transparency=rep(0.5,3), indexAlpha = "all")

Edit:

Only by adjusting the plot.ashape3d function, was I able to remove the edges and vertices:

plot.ashape3d.2 <- function (x, clear = TRUE, col = c(2, 2, 2), byComponents = FALSE,                              indexAlpha = 1, transparency = 1, walpha = FALSE, ...) {  as3d <- x  triangles <- as3d$triang  edges <- as3d$edge  vertex <- as3d$vertex  x <- as3d$x  if (class(indexAlpha) == "character")     if (indexAlpha == "ALL" | indexAlpha == "all")       indexAlpha = 1:length(as3d$alpha)  if (any(indexAlpha > length(as3d$alpha)) | any(indexAlpha <=                                                    0)) {    if (max(indexAlpha) > length(as3d$alpha))       error = max(indexAlpha)    else error = min(indexAlpha)    stop(paste("indexAlpha out of bound : valid range = 1:",                length(as3d$alpha), ", problematic value = ", error,                sep = ""), call. = TRUE)  }  if (clear) {    rgl.clear()  }  if (byComponents) {    components = components_ashape3d(as3d, indexAlpha)    if (length(indexAlpha) == 1)       components = list(components)    indexComponents = 0    for (iAlpha in indexAlpha) {      if (iAlpha != indexAlpha[1])         rgl.open()      if (walpha)         title3d(main = paste("alpha =", as3d$alpha[iAlpha]))      cat("Device ", rgl.cur(), " : alpha = ", as3d$alpha[iAlpha],           "\n")      indexComponents = indexComponents + 1      components[[indexComponents]][components[[indexComponents]] ==                                       -1] = 0      colors = c("#000000", sample(rainbow(max(components[[indexComponents]]))))      tr <- t(triangles[triangles[, 8 + iAlpha] == 2 |                           triangles[, 8 + iAlpha] == 3, c("tr1", "tr2",                                                           "tr3")])      if (length(tr) != 0)         rgl.triangles(x[tr, 1], x[tr, 2], x[tr, 3], col = colors[1 +                                                                    components[[indexComponents]][tr]], alpha = transparency,                       ...)    }  }  else {    for (iAlpha in indexAlpha) {      if (iAlpha != indexAlpha[1])         rgl.open()      if (walpha)         title3d(main = paste("alpha =", as3d$alpha[iAlpha]))      cat("Device ", rgl.cur(), " : alpha = ", as3d$alpha[iAlpha],           "\n")      tr <- t(triangles[triangles[, 8 + iAlpha] == 2 |                           triangles[, 8 + iAlpha] == 3, c("tr1", "tr2",                                                           "tr3")])      if (length(tr) != 0)         rgl.triangles(x[tr, 1], x[tr, 2], x[tr, 3], col = col[1],                       , alpha = transparency, ...)    }  }}alphabunny <- ashape3d(bunny, alpha = c(0.003))plot.ashape3d.2(alphabunny, col=5, indexAlpha = "all", transparency=1)bg3d(1)

The package Rvcg updated to version 0.14 in July 2016, and ball pivoting surface reconstruction was added. The function is vcgBallPivoting:

library(Rvcg) # needs to be >= version 0.14library(rgl)library(onion); data(bunny)# default parametersbunnybp <- vcgBallPivoting(bunny, radius = 0.0022, clustering = 0.2, angle = pi/2)shade3d(bunnybp, col = rainbow(1000), specular = "black")shade3d(bunnybp, col = "pink", specular = "black") # easier to see problem areas.

The ball pivoting and the default parameter settings are not perfect for the Stanford bunny (as noted by cuttlefish44 in the comments radius = 0.0022 does better than the default radius = 0), and you are left with some gaps in the surface. The actual bunny has 2 holes in the base and some scanning limitations contribute to a few other holes (as mentioned here: https://graphics.stanford.edu/software/scanview/models/bunny.html). You may be able to find better parameters, and it's quite fast to use vcgBallPivoting (~0.5 seconds on my machine), but additional effort / methods may be required to close the gaps.