Plotting surface of a sphere using 3d Delaunay triangulated panels in R

Question

[EDIT: more general solutions can be seen in answers to this question]

I'm wondering if anyone can help me plot an approximation of the surface of a sphere using XYZ coordinates. I have tried calculating Delaunay triangulated panels using the package geometry and then to plot iwith rgl. The first attempt, which looks nice, unfortunately created Delaunay 3d triangles that cross through the sphere. I would ultimately like to only plot the surface:

Generate 3d xyz data of a sphere

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)
xyzw <- cbind(x,y,z,w=1)

Calculate 3d Delaunay triangles and plot with rgl:

#install.packages("geometry")
library(geometry)
library(rgl)

tc <- delaunayn(xyzw[,1:3])
open3d()
tetramesh(tc,cbind(x,y,z), alpha=0.2, col=5)
rgl.snapshot("3d_delaunay.png")

Attempt at only returning surface triangles via 2d Delaunay triangulation

tc <- delaunayn(xyzw[,c(1:2)])
open3d()
for(i in seq(nrow(tc))){
    vertices <- c(t(xyzw[tc[i,],]))
    indices <- c( 1, 2, 3)
    shade3d( tmesh3d(vertices,indices) , alpha=0.2, col="cyan")
}
rgl.snapshot("2d_delaunay.png")

Obviously, something is not working. Any help would be greatly appreciated.

Answer 1

What you need is the convex hull, not the Delaunay triangulation. Here is a way with the cxhull package, but you can also use the geometry package.

library(tessellation)

n <- 10
rho <- 1
theta <- seq(0, 2*pi, length.out = n) # azimuthal coordinate running from 0 to 2*pi
phi <- seq(0, pi, length.out = 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)


library(cxhull)
hull <- cxhull(xyz, triangulate = TRUE)

library(rgl)
open3d()
for(facet in hull$facets){
  triangles3d(xyz[facet[["vertices"]], ], color = "green")
}