R语言如何绘制带孔洞的曲面z=f(x,y)？求简便实现方案

需求概述 #

要制作带孔洞的3D图形（用于对比展示原始数据与插值效果），当前通过delaunay包做约束Delaunay三角剖分，再用rgl构建3D网格实现，但想找更高效简便的方法。试过rgl::perps3d.deldir，但该方法会自动填充孔洞，手动移除无效区域的效率极低。另外需要在原始图上叠加黑色的插值区域。

现有数据生成代码 #

x <- y <- seq(-3, 3, by = 0.1)
Grid <- expand.grid(X = x, Y = y)
toremove <- with(Grid, X >= -0.5 & X <= 1 & Y >= -0.5 & Y <= 1)
Grid <- Grid[!toremove, ]
z <- with(Grid, dt(X, df = 1) * dt(Y, df = 1))

两种高效实现方案 #

方案1：rgl + tripack 约束三角剖分 #

tripack的tri.mesh原生支持带边界约束的三角剖分，能直接保留孔洞，不用手动清理无效面，效率更高：

library(rgl)
library(tripack)

# 定义孔洞的闭合边界
hole_border <- rbind(
  c(-0.5, -0.5), c(1, -0.5), c(1, 1), c(-0.5, 1), c(-0.5, -0.5)
)

# 构建带约束的三角网格
tri_mesh <- tri.mesh(Grid$X, Grid$Y, border = hole_border)

# 提取三角面索引并绘制3D网格
open3d()
mesh3d(
  vertices = cbind(Grid$X, Grid$Y, z),
  indices = t(tri_mesh$triangles),
  color = "lightblue"
)

# 生成并叠加黑色插值区域（以线性插值为例）
library(akima)
# 对孔洞区域做插值
interp_data <- interp(
  x = Grid$X, y = Grid$Y, z = z,
  xo = seq(-0.5, 1, by = 0.05),
  yo = seq(-0.5, 1, by = 0.05)
)

# 绘制插值曲面
surface3d(
  interp_data$x, interp_data$y, interp_data$z,
  color = "black", alpha = 0.8
)

方案2：plot3D 蒙版法（最简代码） #

plot3D的surf3D支持通过NA蒙版快速隐藏孔洞区域，完全不用处理三角剖分，代码极简：

library(plot3D)
library(akima)

# 构建完整网格的z矩阵，孔洞区域设为NA
full_grid <- expand.grid(X = x, Y = y)
full_z <- with(full_grid, dt(X, df = 1) * dt(Y, df = 1))
full_z_mat <- matrix(full_z, nrow = length(x), ncol = length(y))
# 蒙版：孔洞位置设为NA
full_z_mat[toremove] <- NA

# 绘制带孔洞的原始曲面
surf3D(x, y, full_z_mat, col = "lightblue", border = "white")

# 生成孔洞区域的插值数据并叠加
interp_hole <- interp(
  Grid$X, Grid$Y, z,
  xo = seq(-0.5, 1, by = 0.05),
  yo = seq(-0.5, 1, by = 0.05)
)

# 添加黑色插值曲面
surf3D(
  interp_hole$x, interp_hole$y, interp_hole$z,
  col = "black", add = TRUE, alpha = 0.8
)

方案对比 #

• 方案1适合需要精细控制三角网格的场景，rgl的交互性更好；
• 方案2代码最简洁，适合快速生成静态可视化，学习成本低；
• 两种方案都能避免perps3d.deldir的孔洞填充问题，效率远超手动清理。

内容的提问来源于stack exchange，提问作者Stéphane Laurent