用xyz坐标绘制三维表面图

我希望有经验的人可以帮助如何从xyz数据准备形状文件。 一个精心准备的数据集的一个很好的例子可以在这里看到彗星Churyumov-Gerasimenko,虽然没有提供创build形状文件的前面的步骤。

我试图更好地理解如何将一个曲面应用到给定的一组XYZ坐标。 使用笛卡尔坐标直接与R包“rgl”,但形状环绕似乎更困难。 我find了R程序包geometry ,它提供了一个QHULL函数的接口。 我试图用这个来计算Delaunay三angular面,然后我可以在rglrgl 。 我无法找出与函数delaunayn相关的一些选项,可能控制这些方面计算的最大距离。 我希望这里的某个人可能会对xyz数据的表面build设有所改进。

使用“斯坦福兔子”数据集的示例:

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

在这里输入图像说明

这个答案让我相信Qhu的R实现可能会有问题。 另外,我现在已经尝试了各种设置(例如delaunayn(bunny, options="Qt") ),效果不大。 在这里列出了 Qhull选项

编辑:

这是一个额外的(更简单的)球体的例子。 即使在这里,面的计算并不总是find最近的相邻顶点(如果你旋转球,你会看到一些穿过内部的方面)。

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

在这里输入图像说明

这是一个使用核密度估计和contour3d函数的misc3d 。 我玩了一段时间,直到我find了一个体面的工作价值。 这并不十分精确,但是你可以调整一些东西来获得更好,更精确的表面。 如果你有超过8GB的内存,那么你可能会增加超出我在这里所做的。

 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) 

在这里输入图像说明在这里输入图像说明

右边的图像是从底部看的,只是为了显示另一个视图。

您可以在kde3d使用较大的n值,但需要更长的时间,并且如果数组变得太大,则可能会用完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") 

在这里输入图像说明


存在更好的,更高效的表面重build方法(例如功率壳,泊松曲面重构,球轴algorithm),但是我不知道R中有任何实现。

这里有一个相关的堆栈溢出文章,提供了一些很好的信息和链接来检出(包括代码链接): 从三维点云表面重build鲁棒algorithm? 。

我想使用alphashape3d软件包find了一个可能的解决scheme。 为了得到一个alpha的可接受值,我必须玩一些,这与给定数据集中的距离有关(例如bunny sd给了我一些见解)。 我仍然试图找出如何更好地控制顶点和边的线条宽度,以免占据主导地位,但这可能与rgl设置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) 

在这里输入图像说明

软件包Rvcg于2016年7月更新至版本0.14,并增加了球形旋转表面重build。 函数是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. 

在这里输入图像说明 在这里输入图像说明

球摆动和默认的参数设置不是完美的斯坦福兔子(如评论radius = 0.0022墨鱼注意到比默认radius = 0 ),并留下一些表面的差距。 实际的兔子在底部有2个孔,一些扫描限制有助于其他几个孔(如这里所述: https : //graphics.stanford.edu/software/scanview/models/bunny.html )。 您可以find更好的参数,使用vcgBallPivoting (在我的机器上约0.5秒)是相当快的,但是可能需要额外的努力/方法来缩小差距。