R语言：颜色随z轴渐变的3D图像

R语言通常被公认的是优秀的绘图能力，但是入坑一年最大的感受是，大部分包着力于2D作图，比如最喜欢的ggplot2，而且这些方面的资料网上也非常丰富；但是3D作图方面，却发现网上资料不是那么多了。

今天要尝试的是绘制一个二元函数的3D图像，重点在于实现图像颜色随z轴渐变。效果如下：

通过persp函数绘制

以上分别通过persp和persp3d两个函数绘制。persp来自R的基础包，而persp3d来自3D绘图包rgl，两者的处理非常相似。

# 输入二元函数,并生成对应的点集
x <- seq(-3, 3, length=50)
y <- seq(-2, 2, length=50)
f <- function(x,y)
{
  return(x^4/2 + y*x^2 - x/2 + y/4)
}
z <- outer(x, y, f)

# expression函数生成标题的函数表达式字符串
title <- expression(f(x,y)==x^4/2 + y*x^2 - x/2 + y/4)

# 通过colorRampPalette函数生成中间色号
nbcol <- 50   #每种颜色间渐变色的数量
jet.colors1 <- colorRampPalette( c("#0000ff", "#00ff00") ) 
color1 <- jet.colors1(nbcol)
jet.colors2 <- colorRampPalette( c("#00ff00", "#ffff00") ) 
color2 <- jet.colors2(nbcol)
jet.colors3 <- colorRampPalette( c("#ffff00", "#cc0000") ) 
color3 <- jet.colors3(nbcol)
color<-matrix(c(color1, color2, color3), ncol=1)

# 计算中心的z值
zfacet <- z[-1, -1] + z[-1, -length(y)] + z[-length(x), -1] + z[-length(x), -length(y)]
# 将z轴值重新编码为颜色索引
facetcol <- cut(zfacet, breaks=nbcol*3)

# persp绘图
par(bg = "white")
persp(x, y, z, col = color[facetcol],
      theta = 45, phi =30, expand = 0.5, 
      xlab = "x", ylab = "y", zlab = "z", zlim = c(-20,60),
      box=TRUE, axes=TRUE,
      ltheta = 120, shade = 0.5, ticktype = "detailed",
      main=title) 
} 


# 以下则是通过persp3d绘图

library("rgl")
# 注意这里的向量处理方式与persp不同
vertcol <- cut(z, nbcol*3)
persp3d(x, y, z,
        col=color[vertcol],smooth=FALSE,lit=FALSE,
        xlab="x", ylab="y", zlab="z", zlim = c(-20,60)
        #,main=title
        ## Herer is a bug with title
        )
# 显示网格线
grid3d(c("x+", "y+", "z"))

颜色处理的过程基本一致，不过中间的细微差别大家应该注意到了，facetcol与vertcol这两个向量分别由zfacet和z矩阵cut来的，而zfacet经过了四角余子式求和的操作，如果在persp中省略了这一步操作你会得到一个蜜汁图像。

具体的原因嘛，可以看上面的第一幅图，注意到surface被分成很多小格吗？绘3D图用的是一个50x50矩阵z的数据，而数据中的2500个点构成了小格的四角，而上色实际是给小格上色的，所以一共49x49个小格，需要将z删掉一行一列，上面余子式相加的操作实际类似于取平均。你当然也可以使用zfacet <- z[-1, -1]试试看，结果没有太大影响，不过显然取平均更严谨一些。然后呢通过cut函数将zfacet划分为若干区间，最后通过color[facetcol]就是按照区间匹配color中的颜色。很神奇吧？你可以用下面的Demo快速演示一下~

Demo: