r - 带有 xyz 坐标的 3d 曲面图

Question

我希望有经验的人可以帮助如何从 xyz 数据准备形状文件。在这里可以看到一个为 Churyumov–Gerasimenko 彗星准备好的数据集的一个很好的例子，尽管没有提供前面创建形状文件的步骤。

我试图更好地理解如何将表面应用于给定的一组 XYZ 坐标。使用 R 包“rgl”直接使用笛卡尔坐标，但是环绕的形状似乎更困难。我找到了 R 包geometry，它提供了QHULL函数的接口。我尝试使用它来计算 Delaunay 三角面，然后我可以在rgl. 我无法弄清楚与该函数相关的一些选项，delaunayn以可能控制计算这些方面的最大距离。我希望这里的某个人可能对从 xyz 数据改进表面构造有一些想法。

使用“Stanford bunnny”数据集的示例：

library(onion)
library(rgl)
library(geometry)
data(bunny)

#XYZ point plot
open3d()
points3d(bunny, col=8, size=0.1)
#rgl.snapshot("3d_bunny_points.png")

#Facets following Delaunay triangulation
tc.bunny <- delaunayn(bunny)
open3d()
tetramesh(tc.bunny, bunny, alpha=0.25, col=8)
#rgl.snapshot("3d_bunny_facets.png")

在此处输入图像描述

这个答案让我相信 Qhull 的 R 实现可能存在问题。另外，我现在尝试了各种设置（例如delaunayn(bunny, options="Qt")），但效果不大。此处概述了 Qhull 选项

编辑：

这是球体的另一个（更简单）示例。即使在这里，面的计算也并不总是找到最近的相邻顶点（如果你旋转球，你会看到一些面穿过内部）。

library(rgl)
library(geometry)
set.seed(1)
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)

set.seed(1)
xyz <- cbind(x,y,z)
tbr = t(surf.tri(xyz, delaunayn(xyz)))
open3d()
rgl.triangles(xyz[tbr,1], xyz[tbr,2], xyz[tbr,3], col = 5, alpha=0.5)
rgl.snapshot("ball.png")

在此处输入图像描述

score 18 · Accepted Answer

这是一种使用核密度估计和contour3d函数的方法misc3d。我一直在玩，直到我找到一个合适的价值levels。它并不完全精确，但您可以调整一些东西以获得更好、更准确的表面。如果你有超过 8GB 的内存，那么你可能会比n我在这里做的更多。

library(rgl)
library(misc3d)
library(onion); data(bunny)

# the larger the n, the longer it takes, the more RAM you need
bunny.dens <- kde3d(bunny[,1],bunny[,2],bunny[,3], n=150, 
    lims=c(-.1,.2,-.1,.2,-.1,.2)) # I chose lim values manually

contour3d(bunny.dens$d, level = 600, 
    color = "pink", color2 = "green", smooth=500)
rgl.viewpoint(zoom=.75)

在此处输入图像描述

右边的图像是从底部开始的，只是为了显示另一个视图。

n您可以为in使用更大的值，kde3d但它会花费更长的时间，并且如果数组变得太大，您可能会耗尽 RAM。您也可以尝试不同的带宽（此处使用默认值）。我从R-Feng & Tierney 2008 中的计算和显示等值面中采用了这种方法。

使用包的非常相似的等值面方法Rvcg：

library(Rvcg)
library(rgl)
library(misc3d)
library(onion); data(bunny)

bunny.dens <- kde3d(bunny[,1],bunny[,2],bunny[,3], n=150, 
    lims=c(-.1,.2,-.1,.2,-.1,.2)) # I chose lim values manually

bunny.mesh <- vcgIsosurface(bunny.dens$d, threshold=600)
shade3d(vcgSmooth(bunny.mesh,"HC",iteration=3),col="pink") # do a little smoothing

在此处输入图像描述

由于它是一种基于密度估计的方法，我们可以通过增加兔子的密度来获得更多收益。我n=400这里也用。成本是计算时间的显着增加，但生成的表面更好：

bunny.dens <- kde3d(rep(bunny[,1], 10), # increase density.
                    rep(bunny[,2], 10),
                    rep(bunny[,3], 10), n=400, 
                    lims=c(-.1,.2,-.1,.2,-.1,.2))

bunny.mesh <- vcgIsosurface(bunny.dens$d, threshold=600)
shade3d(vcgSmooth(bunny.mesh,"HC",iteration=1), col="pink")

在此处输入图像描述

存在更好、更有效的表面重建方法（例如功率外壳、泊松表面重建、球枢轴算法），但我不知道在 R 中已经实现了任何方法。

这是一个相关的 Stack Overflow 帖子，其中包含一些重要信息和要查看的链接（包括代码链接）： robust algorithm for surface rebuild from 3D point cloud? .

score 13 · Accepted Answer

我认为已经找到了一种使用该alphashape3d软件包的可能解决方案。我不得不尝试一下以获得可接受的值alpha，这与给定数据集中的距离有关（例如sd，bunny给了我一些见解）。我仍在试图弄清楚如何更好地控制顶点和边的线宽，以免主导情节，但这可能与rgl.

例子：

library(onion)
library(rgl)
library(geometry)
library(alphashape3d)

data(bunny)
apply(bunny,2,sd)
alphabunny <- ashape3d(bunny, alpha = 0.003)
bg3d(1)
plot.ashape3d(alphabunny, col=c(5,5,5), lwd=0.001, size=0, transparency=rep(0.5,3), indexAlpha = "all")

在此处输入图像描述

编辑：

只有通过调整plot.ashape3d函数，我才能删除边缘和顶点：

plot.ashape3d.2 <- function (x, clear = TRUE, col = c(2, 2, 2), byComponents = FALSE, 
                             indexAlpha = 1, transparency = 1, walpha = FALSE, ...) 
{
  as3d <- x
  triangles <- as3d$triang
  edges <- as3d$edge
  vertex <- as3d$vertex
  x <- as3d$x
  if (class(indexAlpha) == "character") 
    if (indexAlpha == "ALL" | indexAlpha == "all") 
      indexAlpha = 1:length(as3d$alpha)
  if (any(indexAlpha > length(as3d$alpha)) | any(indexAlpha <= 
                                                   0)) {
    if (max(indexAlpha) > length(as3d$alpha)) 
      error = max(indexAlpha)
    else error = min(indexAlpha)
    stop(paste("indexAlpha out of bound : valid range = 1:", 
               length(as3d$alpha), ", problematic value = ", error, 
               sep = ""), call. = TRUE)
  }
  if (clear) {
    rgl.clear()
  }
  if (byComponents) {
    components = components_ashape3d(as3d, indexAlpha)
    if (length(indexAlpha) == 1) 
      components = list(components)
    indexComponents = 0
    for (iAlpha in indexAlpha) {
      if (iAlpha != indexAlpha[1]) 
        rgl.open()
      if (walpha) 
        title3d(main = paste("alpha =", as3d$alpha[iAlpha]))
      cat("Device ", rgl.cur(), " : alpha = ", as3d$alpha[iAlpha], 
          "\n")
      indexComponents = indexComponents + 1
      components[[indexComponents]][components[[indexComponents]] == 
                                      -1] = 0
      colors = c("#000000", sample(rainbow(max(components[[indexComponents]]))))
      tr <- t(triangles[triangles[, 8 + iAlpha] == 2 | 
                          triangles[, 8 + iAlpha] == 3, c("tr1", "tr2", 
                                                          "tr3")])
      if (length(tr) != 0) 
        rgl.triangles(x[tr, 1], x[tr, 2], x[tr, 3], col = colors[1 + 
                                                                   components[[indexComponents]][tr]], alpha = transparency, 
                      ...)
    }
  }
  else {
    for (iAlpha in indexAlpha) {
      if (iAlpha != indexAlpha[1]) 
        rgl.open()
      if (walpha) 
        title3d(main = paste("alpha =", as3d$alpha[iAlpha]))
      cat("Device ", rgl.cur(), " : alpha = ", as3d$alpha[iAlpha], 
          "\n")
      tr <- t(triangles[triangles[, 8 + iAlpha] == 2 | 
                          triangles[, 8 + iAlpha] == 3, c("tr1", "tr2", 
                                                          "tr3")])
      if (length(tr) != 0) 
        rgl.triangles(x[tr, 1], x[tr, 2], x[tr, 3], col = col[1], 
                      , alpha = transparency, ...)
    }
  }
}

alphabunny <- ashape3d(bunny, alpha = c(0.003))
plot.ashape3d.2(alphabunny, col=5, indexAlpha = "all", transparency=1)
bg3d(1)

在此处输入图像描述

score 3 · Accepted Answer

该软件包Rvcg于 2016 年 7 月更新至 0.14 版，并添加了球旋转表面重建。功能是vcgBallPivoting：

library(Rvcg) # needs to be >= version 0.14
library(rgl)
library(onion); data(bunny)

# default parameters
bunnybp <- vcgBallPivoting(bunny, radius = 0.0022, clustering = 0.2, angle = pi/2)
shade3d(bunnybp, col = rainbow(1000), specular = "black")
shade3d(bunnybp, col = "pink", specular = "black") # easier to see problem areas.

球旋转和默认参数设置对于 Stanford bunny 来说并不完美（正如评论中的 cuttlefish44 所指出的radius = 0.0022那样比默认值更好radius = 0），并且您在表面上留下了一些空隙。实际的兔子在底座上有 2 个孔，一些扫描限制会导致其他一些孔（如此处所述：https ://graphics.stanford.edu/software/scanview/models/bunny.html ）。您可能能够找到更好的参数，并且使用起来非常快vcgBallPivoting（在我的机器上大约 0.5 秒），但可能需要额外的努力/方法来缩小差距。

score 0 · Accepted Answer

我目前正在制作一个包RCGAL，它目前提供两种类型的表面重建。它基于RcppCGAL包，该包链接到 C++ 库 CGAL（计算几何算法库）的头文件。

这是斯坦福兔子的高级前表面重建：

这是使用默认参数对斯坦福兔的泊松曲面重建：

这个网格不是很精确。通过使用较小的参数值可以获得更精确的网格spacing：

有关更多信息，请参阅此博客文章。

r - 带有 xyz 坐标的 3d 曲面图

使用“Stanford bunnny”数据集的示例：

编辑：

4 回答 4

例子：

编辑：

Related

Reference