Trisurf Plots in R using Plotly

[This article was first published on R – Modern Data, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

---

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

library(plotly)
library(geometry)

g <- expand.grid(
  u = seq(0, 2 * pi, length.out = 24),
  v = seq(-1, 1, length.out = 8)
)
tp <- 1 + 0.5 * g$v * cos(g$u / 2)
m <- matrix(
  c(tp * cos(g$u), tp * sin(g$u), 0.5 * g$v * sin(g$u / 2)), 
  ncol = 3, dimnames = list(NULL, c("x", "y", "z"))
)

# the key though is running delaunayn on g rather than m
d <- delaunayn(g)
td <- t(d)
#  but using m for plotting rather than the 2d g

# define layout options
axs <- list(
  backgroundcolor="rgb(230,230,230)",
  gridcolor="rgb(255,255,255)",
  showbackground=TRUE,
  zerolinecolor="rgb(255,255,255"
)

# now figure out the colormap
#  start by determining the mean of z for each row
#  of the Delaunay vertices
zmean <- apply(d, MARGIN=1, function(row){mean(m[row,3])})

library(scales)
# result will be slighlty different
#  since colour_ramp uses CIELAB instead of RGB
#  could use colorRamp for exact replication
facecolor = colour_ramp(
  brewer_pal(palette="RdBu")(9)
)(rescale(x=zmean))


plot_ly(
  x = m[, 1], y = m[, 2], z = m[, 3],
  # JavaScript is 0 based index so subtract 1
  i = d[, 1]-1, j = d[, 2]-1, k = d[, 3]-1,
  facecolor = facecolor,
  type = "mesh3d"
) %>%
  layout(
    title="Moebius band triangulation",
    scene=list(xaxis=axs,yaxis=axs,zaxis=axs),
    aspectratio=list(x=1,y=1,z=0.5)
  )

n <- 12
h <- 1/(n-1)
r = seq(h, 1, length.out=n)
theta = seq(0, 2*pi, length.out=36)

g <- expand.grid(r=r, theta=theta)

x <- c(g$r * cos(g$theta),0)
y <- c(g$r * sin(g$theta),0)
z <- sin(x*y)

m <- matrix(
  c(x,y,z), 
  ncol = 3,
  dimnames = list(NULL, c("x", "y", "z"))
)

tri <- delaunayn(m[,1:2])

# now figure out the colormap
zmean <- apply(tri,MARGIN=1,function(row){mean(m[row,3])})

library(scales)
library(rje)
facecolor = colour_ramp(
  cubeHelix(12)
)(rescale(x=zmean))

plot_ly(
  x=x, y=y, z=z,
  i=tri[,1]-1, j=tri[,2]-1, k=tri[,3]-1,
  facecolor=facecolor,
  type="mesh3d"
) %>%
  layout(
    title="Triangulated surface",
    scene=list(
      xaxis=axs,
      yaxis=axs,
      zaxis=axs,
      camera=list(
        eye=list(x=1.75,y=-0.7,z=0.75)
      )
    ),
    aspectratio=list(x=1,y=1,z=0.5)
  )

library(geomorph)
plyFile <- 'http://people.sc.fsu.edu/~jburkardt/data/ply/chopper.ply'
dest <- basename(plyFile)
if (!file.exists(dest)) {
  download.file(plyFile, dest)
}

mesh <- read.ply(dest)
# see getS3method("shade3d", "mesh3d") for details on how to plot 

# plot point cloud
x <- mesh$vb["xpts",]
y <- mesh$vb["ypts",]
z <- mesh$vb["zpts",]
m <- matrix(c(x,y,z), ncol=3, dimnames=list(NULL,c("x","y","z")))

# now figure out the colormap
zmean <- apply(t(mesh$it),MARGIN=1,function(row){mean(m[row,3])})

library(scales)
facecolor = colour_ramp(
  brewer_pal(palette="RdBu")(9)
)(rescale(x=zmean))

plot_ly(
  x = x, y = y, z = z,
  i = mesh$it[1,]-1, j = mesh$it[2,]-1, k = mesh$it[3,]-1,
  facecolor = facecolor,
  type = "mesh3d"
)