studio - superponer graficas en r

Trama de superficie 3D con coordenadas xyz

Question

Espero que alguien con experiencia pueda ayudarme a preparar los archivos shape desde datos xyz. Un gran ejemplo de un conjunto de datos bien preparado se puede ver here para el cometa Churyumov-Gerasimenko, aunque los pasos anteriores en la creación del archivo de forma no se proporcionan.

Estoy intentando comprender mejor cómo aplicar una superficie a un conjunto dado de coordenadas XYZ. El uso de coordenadas cartesianas es sencillo con el paquete R "rgl", pero las formas que se envuelven parecen más difíciles. Encontré la geometry paquete R, que proporciona una interfaz para QHULL funciones QHULL . Intenté usar esto para calcular las facetas triangulares de Delaunay, que luego puedo trazar en rgl . No puedo descubrir algunas de las opciones asociadas con la función delaunayn para posiblemente controlar las distancias máximas que se calculan estas facetas. Espero que alguien aquí tenga algunas ideas para mejorar la construcción de la superficie a partir de datos xyz.

Ejemplo usando el conjunto de datos "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")

Esta respuesta me hace creer que podría haber un problema con la implementación R de Qhull. Además, ahora he probado varias configuraciones (por ejemplo, delaunayn(bunny, options="Qt") ) con poco efecto. Las opciones de Qhull se describen here

Editar:

Aquí hay un ejemplo adicional (más simple) de una esfera. Incluso aquí, el cálculo de las facetas no siempre encuentra los vértices vecinos más cercanos (si giras la bola verás algunas facetas cruzando el 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")

Answer 1

Aquí hay un enfoque que usa la estimación de la densidad del kernel y la función misc3d de misc3d . Jugué hasta que encontré un valor para los levels que funcionaban decentemente. No es perfectamente preciso, pero es posible que pueda modificar las cosas para obtener una superficie mejor y más precisa. Si tiene más de 8 GB de memoria, entonces puede aumentar n más allá de lo que hice aquí.

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)

La imagen de la derecha es desde abajo, solo para mostrar otra vista.

Puede usar un valor mayor para n en kde3d pero llevará más tiempo, y puede quedarse sin RAM si la matriz se vuelve demasiado grande. También podría intentar con un ancho de banda diferente (por defecto se usa aquí). Tomé este enfoque de Computing and Displaying Isosurfaces en R - Feng & Tierney 2008 .

Enfoque isosuperficie muy similar utilizando el paquete 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

Dado que se trata de un enfoque basado en la estimación de la densidad, podemos obtener un poco más de él al aumentar la densidad del conejito. También uso n=400 aquí. El costo es un aumento significativo en el tiempo de cálculo, pero la superficie resultante es una liebre mejor:

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

Existen mejores métodos de reconstrucción de superficie más eficientes (p. Ej., Corteza de poder, reconstrucción de superficie de Poisson, algoritmo de pivote de bola), pero todavía no sé si se han implementado en R.

Aquí hay una publicación relevante de Desbordamiento de pila con información importante y enlaces para revisar (incluidos los enlaces al código): algoritmo robusto para la reconstrucción de la superficie a partir de la nube de puntos 3D? .

Answer 2

Creo que he encontrado una posible solución usando el paquete alphashape3d . Tuve que jugar un poco para obtener un valor aceptable para alpha , que está relacionado con las distancias en el conjunto de datos dado (por ejemplo, sd de bunny me dio una idea). Todavía intento descubrir cómo controlar mejor el ancho de las líneas en vértices y bordes para no dominar la trama, pero esto probablemente esté relacionado con la configuración en rgl .

Ejemplo:

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

Editar:

Solo al ajustar la función plot.ashape3d , pude eliminar los bordes y vértices:

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)

Answer 3

El paquete Rvcg actualizó a la versión 0.14 en julio de 2016 y se agregó la reconstrucción de la superficie pivotante de la bola. La función es 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.

La rotación de la bola y la configuración predeterminada de los parámetros no son perfectos para el conejito de Stanford (como lo notó la sepia44 en el radius = 0.0022 comentarios radius = 0.0022 es mejor que el radius = 0 predeterminado radius = 0 ), y le quedan algunas lagunas en la superficie. El conejito real tiene 2 agujeros en la base y algunas limitaciones de escaneo contribuyen a algunos otros agujeros (como se menciona aquí: https://graphics.stanford.edu/software/scanview/models/bunny.html ). Es posible que pueda encontrar mejores parámetros, y es bastante rápido usar vcgBallPivoting (~ 0.5 segundos en mi máquina), pero es posible que se requieran esfuerzos / métodos adicionales para cerrar los espacios.