Трехмерная поверхность с координатами xyz

Я надеюсь, что кто-то с опытом может помочь в том, как вы готовите файлы формы из xyz-данных. Хороший пример хорошо подготовленного набора данных можно увидеть здесь здесь для кометы Чурюмов-Герасименко, хотя предыдущие шаги по созданию файл формы не предоставляются.

Я пытаюсь лучше понять, как применить поверхность к заданному набору координат XYZ. Использование декартовых координат прямо с R-пакетом "rgl", однако формы, которые обертываются, кажутся более трудными. Я нашел пакет R geometry, который предоставляет интерфейс для функций QHULL. Я попытался использовать это для вычисления треугольных граней Делоне, которые затем можно построить в rgl. Я не могу определить некоторые параметры, связанные с функцией delaunayn, чтобы, возможно, управлять максимальными расстояниями, которые вычисляются этими гранями. Я надеюсь, что у кого-то здесь могут быть некоторые идеи по улучшению построения поверхности из данных xyz.

Пример с использованием набора данных "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")

Этот ответ заставляет меня поверить, что может возникнуть проблема с реализацией R Qhull. Кроме того, я сейчас пробовал различные настройки (например, delaunayn(bunny, options="Qt")) с небольшим эффектом. Опции Qhull выделены здесь

Edit:

Вот еще один (более простой) пример сферы. Даже здесь вычисление фасетов не всегда находит ближайшие соседние вершины (если вы поворачиваете шар, вы увидите некоторые грани, пересекающие внутренность).

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

Ответ 1

Ниже приведен подход с использованием оценки плотности ядра и функции contour3d от misc3d. Я играл, пока не нашел значение для levels, которое работало прилично. Это не совсем точно, но вы можете настроить все, чтобы получить лучшую, более точную поверхность. Если у вас более 8 ГБ памяти, вы можете увеличить 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)

Изображение справа снизу, просто чтобы показать другое представление.

Вы можете использовать большее значение для n в kde3d, но это займет больше времени, и может закончиться нехватка ОЗУ, если массив становится слишком большим. Вы также можете попробовать другую полосу пропускания (по умолчанию используется здесь). Я использовал этот подход из Вычисление и отображение изоповерхностей в 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")

Лучше, существуют более эффективные методы реконструкции поверхности (например, силовая коса, реконструкция поверхности Пуассона, алгоритм шарового шарнира), но я не знаю, что какие-либо были реализованы в R.

Здесь соответствующая запись с отличной информацией и ссылками для проверки (включая ссылки на код): надежный алгоритм восстановления поверхности из облака 3D-точек?.

Ответ 2

Я думаю, нашел одно возможное решение, используя пакет alphashape3d. Мне пришлось немного поиграть, чтобы получить приемлемое значение для alpha, которое связано с расстояниями в данном наборе данных (например, sd of bunny). Я все еще пытаюсь понять, как лучше контролировать ширину линий в вершинах и краях, чтобы не доминировать над сюжетом, но это, вероятно, связано с настройками в 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")

Edit:

Только путем настройки функции 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)

Ответ 3

Пакет Rvcg обновлен до версии 0.14 в июле 2016 года и добавлена реконструкция поверхности шарового шарнира. Функция 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.

Шаблон поворота и настройки параметров по умолчанию не идеальны для кролика Стэнфорда (как отмечено cuttlefish44 в комментариях radius = 0.0022 лучше, чем по умолчанию radius = 0), и вы остаетесь с некоторыми пробелами на поверхности. Фактический кролик имеет 2 отверстия в основании, а некоторые ограничения сканирования вносят вклад в несколько других отверстий (как указано здесь: https://graphics.stanford.edu/software/scanview/models/bunny.html). Вы можете найти лучшие параметры, и довольно быстро использовать vcgBallPivoting (~ 0,5 секунды на моей машине), но для закрытия пробелов может потребоваться дополнительное усилие/методы.