与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")

Answer 1

这是一种使用核密度估计和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,kde3d但是它需要更长的时间,如果数组太大,你可能会用完RAM.您还可以尝试不同的带宽(此处默认使用).我在R-Feng&Tierney 2008中采用了计算和显示等值面的方法.

使用Rvcg包的非常相似的isosurface方法:

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帖子,其中包含一些很棒的信息和链接(包括代码链接): 用于从3D点云进行曲面重建的强大算法？.

Answer 2

我认为使用该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)