用xyz坐标绘制三维表面图
我希望有经验的人可以帮助如何从xyz数据准备形状文件。 一个精心准备的数据集的一个很好的例子可以在这里看到彗星Churyumov-Gerasimenko，虽然没有提供创build形状文件的前面的步骤。
我试图更好地理解如何将一个曲面应用到给定的一组XYZ坐标。 使用笛卡尔坐标直接与R包“rgl”，但形状环绕似乎更困难。 我find了R程序包geometry ，它提供了一个QHULL函数的接口。 我试图用这个来计算Delaunay三angular面，然后我可以在rgl中rgl 。 我无法找出与函数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秒）是相当快的，但是可能需要额外的努力/方法来缩小差距。
