I am hoping someone with experience can help in how one prepares the shape files from xyz data. A great example of a well-prepared dataset can be seen here for the comet Churyumov–Gerasimenko, although the preceding steps in creating the shape file are not provided.
I'm trying to better understand how to apply a surface to a given set of XYZ coordinates. Using Cartesian coordinates is straight forward with the R package "rgl", however shapes that wrap around seem more difficult. I found the R package
#XYZ point plot
points3d(bunny, col=8, size=0.1)
#Facets following Delaunay triangulation
tc.bunny <- delaunayn(bunny)
tetramesh(tc.bunny, bunny, alpha=0.25, col=8)
n <- 10
rho <- 1
theta <- seq(0, 2*pi,, n) # azimuthal coordinate running from 0 to 2*pi
phi <- seq(0, pi,, n) # polar coordinate running from 0 to pi (colatitude)
grd <- expand.grid(theta=theta, phi=phi)
x <- rho * cos(grd$theta) * sin(grd$phi)
y <- rho * sin(grd$theta) * sin(grd$phi)
z <- rho * cos(grd$phi)
xyz <- cbind(x,y,z)
tbr = t(surf.tri(xyz, delaunayn(xyz)))
rgl.triangles(xyz[tbr,1], xyz[tbr,2], xyz[tbr,3], col = 5, alpha=0.5)
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 need bunny.dens <- kde3d(bunny[,1],bunny[,2],bunny[,3], n=150, lims=c(-.1,.2,-.1,.2,-.1,.2)) # I chose lim values manually contour3d(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
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
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 manually bunny.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?.