Tracé de surface 3D avec des coordonnées xyz

J’espère que quelqu’un avec de l’expérience peut vous aider à préparer les fichiers de forme à partir de données xyz. Un bon exemple de jeu de données bien préparé peut être vu ici pour la comète Churyumov – Gerasimenko, bien que les étapes précédentes de création du fichier de forme ne soient pas fournies.

J’essaie de mieux comprendre comment appliquer une surface à un ensemble donné de coordonnées XYZ. L’utilisation des coordonnées cartésiennes est simple avec le package R “rgl”, cependant les formes qui semblent envelopper semblent plus difficiles. J’ai trouvé la geometry paquet R, qui fournit une interface aux fonctions QHULL . J’ai essayé de l’utiliser pour calculer les facettes sortingangulées de Delaunay, que je peux ensuite tracer en rgl . Je n’arrive pas à comprendre certaines des options associées à la fonction delaunayn pour contrôler éventuellement les distances maximales calculées par ces facettes. J’espère que quelqu’un ici pourrait avoir des idées pour améliorer la construction de surface à partir de données xyz.

Exemple d’utilisation du jeu de données “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 sortingangulation tc.bunny <- delaunayn(bunny) open3d() tetramesh(tc.bunny, bunny, alpha=0.25, col=8) #rgl.snapshot("3d_bunny_facets.png") 

entrer la description de l'image ici

Cette réponse me fait croire qu’il pourrait y avoir un problème avec la mise en œuvre de R de Qhull. En outre, j’ai maintenant essayé divers parameters (par exemple, delaunayn(bunny, options="Qt") ) avec peu d’effet. Les options Qhull sont décrites ici

Modifier:

Voici un exemple supplémentaire (plus simple) d’une sphère. Même ici, le calcul des facettes ne trouve pas toujours les sumts voisins les plus proches (si vous faites pivoter la balle, vous verrez certaines facettes traverser l’intérieur).

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

entrer la description de l'image ici

Voici une approche utilisant l’estimation de la densité du kernel et la fonction misc3d de misc3d . J’ai joué jusqu’à ce que je trouve une valeur pour des levels fonctionnant convenablement. Ce n’est pas parfaitement précis, mais vous pourrez peut-être modifier les choses pour obtenir une surface meilleure et plus précise. Si vous avez plus de 8 Go de mémoire, vous pourrez peut-être augmenter ce que j’ai fait ici.

 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) 

entrer la description de l'image icientrer la description de l'image ici

L'image à droite est en bas, juste pour montrer une autre vue.

Vous pouvez utiliser une valeur plus grande pour n dans kde3d mais cela prendra plus de temps et vous risquez de manquer de mémoire vive si la baie devient trop volumineuse. Vous pouvez également essayer une bande passante différente (utilisée par défaut ici). J'ai pris cette approche de l ' informatique et de l' affichage des isosurfaces dans R - Feng & Tierney 2008 .


Approche isosurface très similaire utilisant le 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 

entrer la description de l'image ici

Puisqu'il s'agit d'une approche basée sur l'estimation de la densité, nous pouvons en tirer un peu plus en augmentant la densité du lapin. J'utilise aussi n=400 ici. Le coût est une augmentation significative du temps de calcul, mais la surface résultante est mieux:

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

entrer la description de l'image ici


Il existe des méthodes de reconstruction de surface meilleures et plus efficaces (par exemple, croûte de puissance, reconstruction de surface de Poisson, algorithme de pivot de boule), mais je ne sais pas si elles ont été implémentées dans R, pour le moment.

Voici un article pertinent sur Stack Overflow avec des informations et des liens intéressants à extraire (y compris des liens vers du code): algorithme robuste pour la reconstruction de surface à partir d'un nuage de points 3D? .

Je pense avoir trouvé une solution possible en utilisant le package alphashape3d . J’ai dû jouer un peu pour obtenir une valeur acceptable pour l’ alpha , qui est liée aux distances dans l’dataset donné (par exemple, sd de bunny m’a donné un aperçu). J’essaie toujours de comprendre comment mieux contrôler la largeur des lignes dans les sumts et les arêtes afin de ne pas dominer le tracé, mais cela est probablement lié aux parameters de rgl .

Exemple:

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

entrer la description de l'image ici

Modifier:

Ce n'est qu'en ajustant la fonction plot.ashape3d j'ai pu supprimer les arêtes et les sumts:

 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) 

entrer la description de l'image ici

Le package Rvcg mis à jour vers la version 0.14 en juillet 2016 et la reconstruction de la surface de pivotement des billes ont été ajoutés. La fonction est 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. 

entrer la description de l'image ici entrer la description de l'image ici

Le pivotement de la balle et les parameters par défaut ne sont pas parfaits pour le lapin de Stanford (comme le signale radius = 0.0022 dans le radius = 0.0022 commentaires radius = 0.0022 fait mieux que le radius = 0 par défaut radius = 0 ), il vous rest quelques lacunes dans la surface. Le lapin a deux trous dans la base et certaines limitations de balayage consortingbuent à quelques autres trous (comme mentionné ici: https://graphics.stanford.edu/software/scanview/models/bunny.html ). Vous pourrez peut-être trouver de meilleurs parameters, et il est assez rapide d'utiliser vcgBallPivoting (~ 0,5 seconde sur ma machine), mais des efforts / méthodes supplémentaires peuvent être nécessaires pour combler les lacunes.