具有xyz坐标的三维曲面图

I am hoping someone with experience can help in how one prepares the shape files from xyz data. A great example of a well-prepared dataset can be seen here for the comet Churyumov–Gerasimenko, although the preceding steps in creating the shape file are not provided.

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

I'm trying to better understand how to apply a surface to a given set of XYZ coordinates. Using Cartesian coordinates is straight forward with the R package "rgl", however shapes that wrap around seem more difficult. I found the R package geometry, which provides an interface to QHULL functions. I tried using this to calculate Delaunay triangulated facets, which I can then plot in rgl. I'm unable to figure out some of the options associated with the function delaunayn to possibly control the maximum distances that these facets are calculated. I am hoping that someone here might have some ideas on improving the surface construction from xyz data.

我想更好地理解如何将曲面应用到给定的XYZ坐标集合中。使用笛卡尔坐标与R包“rgl”是直接向前的，但是环绕的形状看起来更困难。我找到了R包几何，它提供了QHULL函数的接口。我尝试使用它来计算Delaunay三角剖分，然后可以在rgl中绘图。我无法找出与函数delaunayn相关的一些选项来控制计算这些facet的最大距离。我希望在座的各位能够对改进xyz数据的表面结构有一些想法。

Example using "Stanford bunnny" dataset:

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

This answer makes me believe that there might be a problem with the R implementation of Qhull. Also, I have now tried various settings (e.g. delaunayn(bunny, options="Qt")) with little effect. Qhull options are outlined here

这个答案让我相信Qhull的R实现可能存在问题。此外，我还尝试了各种设置(例如delaunayn(bunny, options=“Qt”)，效果都不太好。Qhull选项在这里概述

Edit:

Here is an additional (simpler) example of a sphere. Even here, the calculation of facets does not always find the closest neighboring vertices (if you rotate the ball you will see some facets crossing through the interior).

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

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

3 个解决方案

#1

Here is an approach using kernel density estimation and the contour3d function from misc3d. I played around until I found a value for levels that worked decently. It's not perfectly precise, but you may be able to tweak things to get a better, more accurate surface. If you have more than 8GB of memory, then you may be able to increase n beyond what I did here.

这里有一种方法使用内核密度估计和轮廓3d函数从misc3d。我一直在玩，直到我找到了一种有效的水平。它不是很精确，但是你可以对它进行调整以得到更好、更精确的表面。如果你有超过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)

The image on the right is from the bottom, just to show another view.

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

You can use a larger value for n in kde3d but it will take longer, and you may run out of RAM if the array becomes too large. You could also try a different bandwidth (default used here). I took this approach from Computing and Displaying Isosurfaces in R - Feng & Tierney 2008.

在kde3d中可以对n使用较大的值，但这需要更长的时间，如果数组太大，可能会耗尽RAM。您也可以尝试不同的带宽(默认使用这里)。我在2008年的R - Feng & Tierney的计算和显示等表面上采取了这种方法。

---

Very similar isosurface approach using the Rvcg package:

非常相似的等表面方法使用的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

Since it's a density estimation based approach, we can get a little more out of it by increasing the density of the bunny. I also use n=400 here. The cost is a significant increase in computation time, but the resulting surface is a hare better:

因为这是一种基于密度估计的方法，我们可以通过增加兔子的密度得到更多。这里也用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")

---

Better, more efficient surface reconstruction methods exist (e.g. power crust, Poisson surface reconstruction, ball-pivot algorithm), but I don't know that any have been implemented in R, yet.

更好，更有效的表面重建方法存在(例如:power crust, Poisson surface reconstruction, ball-pivot算法)，但是我不知道在R中有没有实现过。

Here's a relevant Stack Overflow post with some great information and links to check out (including links to code): robust algorithm for surface reconstruction from 3D point cloud?.

这里有一个相关的栈溢出贴子，包含一些很好的信息和链接(包括到代码的链接):从3D点云进行表面重构的鲁棒算法?

#2

I think have found one possible solution using the alphashape3d package. I had to play around a bit to get an acceptable value for alpha, which is related to the distances in the given data set (e.g. sd of bunny gave me some insight). I'm still trying to figure out how to better control the width of lines in vertices and edges so as not to dominate the plot, but this is probably related to settings in rgl.

我认为已经找到了一种使用alphashape3d包的解决方案。为了得到alpha的可接受值，我不得不做了一点修改，这个值与给定数据集中的距离有关(例如，兔子的sd给了我一些见解)。我仍在试图弄清楚如何更好地控制顶点和边中的线的宽度，以免在图中占主导地位，但这可能与rgl中的设置有关。

Example:

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

Edit:

Only by adjusting the plot.ashape3d function, was I able to remove the edges and vertices:

只有通过调整情节。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)

#3

The package Rvcg updated to version 0.14 in July 2016, and ball pivoting surface reconstruction was added. The function is vcgBallPivoting:

在2016年7月，将Rvcg软件包更新为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.

The ball pivoting and the default parameter settings are not perfect for the Stanford bunny (as noted by cuttlefish44 in the comments radius = 0.0022 does better than the default radius = 0), and you are left with some gaps in the surface. The actual bunny has 2 holes in the base and some scanning limitations contribute to a few other holes (as mentioned here: https://graphics.stanford.edu/software/scanview/models/bunny.html). You may be able to find better parameters, and it's quite fast to use vcgBallPivoting (~0.5 seconds on my machine), but additional effort / methods may be required to close the gaps.

球的旋转和默认参数设置并不适合Stanford bunny(正如在comments radius = 0.0022中的cuttlefish44所指出的那样，比默认的radius = 0要好)，并且会在表面留下一些空隙。实际的兔子在底部有两个洞，一些扫描限制导致了一些其他的洞(如这里提到的:https://graphics.stanford.edu/software/scanview/models/bunny.html)。您可能可以找到更好的参数，并且使用vcgBallPivoting(在我的机器上大约是0.5秒)非常快，但是可能需要额外的努力和方法才能消除差距。