在R中绘制3D等值面失败：如何实现MATLAB patch式3D绘图？

问题背景 #

需要在R中绘制和MATLAB patch(isosurface(...)) 等价的3D形状，但使用plotly的mesh3d时只得到空坐标系，无形状显示。以下是现有R代码和参考的MATLAB代码。

现有R代码（无法正常出图） #

library(pracma)
n1=60
n2=60
n3=60
data = expand.grid(x=linspace(-3,3,n1),y=linspace(-3,3,n2),z=linspace(-3,3,n3))

someData <- rep(0, 60*60*60)  
V <- array(someData, c(60, 60, 60))
count = 1
lambda = c(2,1,1)
for (ii in 1:n1){
    for (jj in 1:n2){
        for (kk in 1:n3){
            V[ii, jj, kk] = t(as.matrix(lambda))%*%as.matrix(c(data$x[count],data$y[count],data$z[count]))
            count = count+1
        }
    }
}

library(reshape2)
library(plotly)
M = melt(V)
plot_ly(type = "mesh3d",
        x = data$x,
        y = data$y,
        z = data$z,
        i = M$Var1,
        j = M$Var2,
        k = M$Var3,
        facecolor = rep(toRGB(colorRampPalette(c("navy", "blue"))(6)), each = 2)
)

参考MATLAB代码 #

n1 = 60;
n2 = 60;
n3 = 60;
[y,x,z] = ndgrid(linspace(-3,3,n1),linspace(-3,3,n2),linspace(-3,3,n3));
V = zeros(n1, n2, n3);
lambda = [2;1;1];
for ii = 1: n1
    for jj = 1:n2
        for kk =1:n3
            V(ii, jj, kk) =  lambda'*[abs(y(ii, jj, kk)); abs(x(ii, jj, kk)); abs(z(ii, jj, kk))];
        end
    end
end
p = patch(isosurface(x,y,z,V,1));
p.FaceColor = 'cyan';
p.EdgeColor = 'none';
view(3);
camlight
axis equal
lighting gouraud
box on

问题分析 #

现有R代码的核心错误：

1. 计算V时缺失abs()：MATLAB代码中对坐标取了绝对值，R代码没加，导致V的取值范围不符合等值面要求。
2. 坐标顺序不匹配：expand.grid的坐标排列和MATLAB的ndgrid不一致，导致V的计算对应关系错误。
3. mesh3d参数错误：i/j/k需要的是等值面的三角面顶点索引，不是原始网格的维度索引，直接用melt(V)的结果完全不对。

解决方案 #

使用misc3d包提取等值面的顶点和三角面信息，再用plotly绘制，步骤如下：

1. 安装依赖包 #

install.packages(c("pracma", "misc3d", "plotly", "reshape2"))

2. 修正数据生成与V的计算 #

library(pracma)
library(misc3d)
library(plotly)

n1=60
n2=60
n3=60
# 用和MATLAB ndgrid一致的坐标生成方式
x <- linspace(-3,3,n1)
y <- linspace(-3,3,n2)
z <- linspace(-3,3,n3)

# 构建V矩阵，加入abs()和MATLAB保持一致
V <- array(0, c(n1, n2, n3))
lambda = c(2,1,1)
for (ii in 1:n1){
    for (jj in 1:n2){
        for (kk in 1:n3){
            V[ii, jj, kk] = lambda %*% abs(c(x[ii], y[jj], z[kk]))
        }
    }
}

3. 提取等值面数据 #

# 提取等值为1的面（和MATLAB一致），draw=FALSE表示不直接绘图
iso_data <- contour3d(V, level = 1, x = x, y = y, z = z, draw = FALSE)

4. 用plotly绘制3D形状 #

# plotly的mesh3d索引从0开始，R的索引从1开始，所以要减1
vertices <- iso_data$v
faces <- iso_data$triangles - 1

plot_ly(type = "mesh3d",
        x = vertices[,1],
        y = vertices[,2],
        z = vertices[,3],
        i = faces[,1],
        j = faces[,2],
        k = faces[,3],
        facecolor = toRGB("cyan"), # 和MATLAB的FaceColor一致
        edgecolor = toRGB("none"), # 隐藏边线
        opacity = 1) %>%
  layout(scene = list(aspectmode = "data", # 对应MATLAB的axis equal
                      xaxis = list(title = "x"),
                      yaxis = list(title = "y"),
                      zaxis = list(title = "z")))

效果说明 #

运行上述代码后，会得到和MATLAB代码一致的3D形状：一个沿x轴拉伸的菱形多面体，面为青色，无边线，坐标系比例均等。

内容的提问来源于stack exchange，提问作者Sparsity