R + построено: твердое тело вращения

У меня есть функция r(x), что я хочу вращаться вокруг x оси, чтобы получить тело вращения, что я хотел бы добавить к существующему plot_ly участка с использованием add_surface (окрасили x).

Вот пример:

library(dplyr)
library(plotly)

# radius depends on x
r <- function(x) x^2

# interval of interest
int <- c(1, 3)

# number of points along the x-axis
nx <- 20

# number of points along the rotation
ntheta <- 36

# set x points and get corresponding radii
coords <- data_frame(x = seq(int[1], int[2], length.out = nx), r = r(x))

# for each x: rotate r to get y and z coordinates
# edit: ensure 0 and pi are both amongst the angles used
coords %<>%
  rowwise() %>%
  do(data_frame(x = .$x, r = .$r,
                theta = seq(0, pi, length.out = ntheta / 2 + 1) %>%
                c(pi + .[-c(1, length(.))]))) %>%

  ungroup %>%
  mutate(y = r * cos(theta), z = r * sin(theta))

# plot points to make sure the coordinates define the desired shape
coords %>%
  plot_ly(x = ~x, y = ~y, z = ~z, color = ~x) %>%
  add_markers()

3D scatter plot

Как я могу сформировать фигуру, указанную выше, как поверхность plotly (идеально открытую с обоих концов)?


Редактировать (1):

Вот моя лучшая попытка:

# get all x & y values used (sort to connect halves on the side)
xs <-
  unique(coords$x) %>%
  sort
ys <-
  unique(coords$y) %>%
  sort

# for each possible x/y pair: get z^2 value
coords <-
  expand.grid(x = xs, y = ys) %>%
  as_data_frame %>%
  mutate(r = r(x), z2 = r^2 - y^2)

# format z coordinates above x/y plane as matrix where columns
# represent x and rows y
zs <- matrix(sqrt(coords$z2), ncol = length(xs), byrow = TRUE)

# format x coordiantes as matrix as above (for color gradient)
gradient <-
  rep(xs, length(ys)) %>%
  matrix(ncol = length(xs), byrow = TRUE)

# plot upper half of shape as surface
p <- plot_ly(x = xs, y = ys, z = zs, surfacecolor = gradient,
             type = "surface", colorbar = list(title = 'x'))

# plot lower have of shape as second surface
p %>%
  add_surface(z = -zs, showscale = FALSE)

two 3D surfaces attempt

Хотя это дает желаемую форму,

  1. Он имеет "зубы бритвы" близко к плоскости x/y.
  2. Части половины не касаются. (разрешен путем включения 0 и pi в theta векторы)
  3. Я не понял, как покрасить его вместо x вместо z (хотя я пока не очень разбираюсь в этом). (разрешено gradient матрицей)

edit (2):

Вот попытка использовать одну поверхность:

# close circle in y-direction
ys <- c(ys, rev(ys), ys[1])

# get corresponding z-values
zs <- rbind(zs, -zs[nrow(zs):1, ], zs[1, ])

# as above, but for color gradient
gradient <-
  rbind(gradient, gradient[nrow(gradient):1, ], gradient[1, ])

# plot single surface
plot_ly(x = xs, y = ys, z = zs, surfacecolor = gradient,
        type = "surface", colorbar = list(title = 'x'))

Удивительно, но в то время как это должно соединить две половины, ортогональные плоскости x/y две создают полную форму, она по-прежнему страдает от того же эффекта "бритвенных зубов", что и выше:

single 3D surface attempt


edit (3):

Оказывается, недостающие части возникают из-за того, что z -values является NaN когда он близок к 0:

# color points 'outside' the solid purple
gradient[is.nan(zs)] <- -1

# show those previously hidden points
zs[is.nan(zs)] <- 0

# plot exactly as before
plot_ly(x = xs, y = ys, z = zs, surfacecolor = gradient,
        type = "surface", colorbar = list(title = 'x'))

NaN-override attempt

Это может быть вызвано численной неустойчивостью подстановки, когда r^2 и y слишком близки, что приводит к отрицательному вводу для sqrt где фактический ввод все еще неотрицателен.

Эти швы не связаны с численными проблемами, так как даже при рассмотрении +-4 "близко" к нулю эффект "бритвенных зубов" нельзя полностью исключить:

# re-calculate z-values rounding to zero if 'close'
eps <- 4
zs <- with(coords, ifelse(abs(z2) < eps, 0, sqrt(z2))) %>%
      matrix(ncol = length(xs), byrow = TRUE) %>%
      rbind(-.[nrow(.):1, ], .[1, ])

# plot exactly as before
plot_ly(x = xs, y = ys, z = zs, surfacecolor = gradient,
        type = "surface", colorbar = list(title = 'x'))

eps-attempt

Ответы

Ответ 1

Одним из решений было бы перевернуть ваши оси так, чтобы вы вращались вокруг оси z, а не оси x. Я не знаю, возможно ли это, учитывая существующую диаграмму, на которую вы добавляете эту цифру, но она легко решает проблему "зубов".

xs <- seq(-9,9,length.out = 20)
ys <- seq(-9,9,length.out = 20)

coords <-
  expand.grid(x = xs, y = ys) %>%
  mutate(z2 = (x^2 + y^2)^(1/4))

zs <- matrix(coords$z2, ncol = length(xs), byrow = TRUE)

plot_ly(x = xs, y = ys, z = zs, surfacecolor = zs,
             type = "surface", colorbar = list(title = 'x')) %>% 
  layout(scene = list(zaxis = list(range = c(1,3))))

enter image description here

Ответ 2

интересный вопрос, я изо всех сил пытался использовать поверхностную плотность для улучшения вашего решения. Существует хак, который вы можете сделать с разбиением на несколько строк, что может показаться приятным для этого. Например, только изменения, сделанные в оригинале, включают использование большего количества x точек: от nx до 1000 и изменение add_markers в add_lines. Не может быть масштабируемым, но отлично работает для такого размера данных :)

library(dplyr)
library(plotly)

# radius depends on x
r <- function(x) x^2

# interval of interest
int <- c(1, 3)

# number of points along the x-axis
nx <- 1000

# number of points along the rotation
ntheta <- 36

# set x points and get corresponding radii
coords <- data_frame(x = seq(int[1], int[2], length.out = nx), r = r(x))

# for each x: rotate r to get y and z coordinates
# edit: ensure 0 and pi are both amongst the angles used
coords %<>%
  rowwise() %>%
  do(data_frame(x = .$x, r = .$r,
                theta = seq(0, pi, length.out = ntheta / 2 + 1) %>%
                  c(pi + .[-c(1, length(.))]))) %>%

  ungroup %>%
  mutate(y = r * cos(theta), z = r * sin(theta))

# plot points to make sure the coordinates define the desired shape
coords %>%
  plot_ly(x = ~x, y = ~y, z = ~z, color = ~x) %>%
  add_lines()

enter image description here

Лучший, Джонни

Ответ 3

У меня была еще одна трещина и у меня было более близкое решение, используя "поверхностный" тип. Что помогло посмотреть результаты вашего первого участка поверхности с nx = 5 и ntheta = 18. Причина, по которой это извращение, объясняется тем, как она связывает столбцы в zs (через точки x). Он должен связываться с частичным движением вокруг большего кольца вокруг него, и это заставляет плотность всплывать, чтобы соответствовать этой точке.

Я не могу избавиться от этого неуверенного поведения 100%. Я внес эти изменения:

  1. добавьте несколько небольших точек в тету по краям: где две плотности объединены. Это уменьшает размер неустойчивой части, так как есть еще несколько точек, близких к границе
  2. вычисление до mod zs до zs2: убедитесь, что каждое кольцо имеет равную размерность с внешним кольцом, добавив 0 in.
  3. увеличение nx до 40 и уменьшение ntheta до 18 - больше x делает шаг меньше. уменьшить ntheta для времени выполнения, как я добавил на большее количество очков

шаги приходят в том, как он пытается объединить х-кольца. Теоретически, если у вас больше х колец, это должно устранить эту неровность, но это займет много времени.

Я не думаю, что это отвечает на Q 100%, и я не уверен, что эта библиотека лучше всего подходит для этой работы. Свяжитесь с нами, если у вас есть Q.

library(dplyr)
library(plotly)

# radius depends on x
r <- function(x) x^2

# interval of interest
int <- c(1, 3)

# number of points along the x-axis
nx <- 40

# number of points along the rotation
ntheta <- 18

# set x points and get corresponding radii
coords <- data_frame(x = seq(int[1], int[2], length.out = nx), r = r(x))

# theta: add small increments at the extremities for the density plot
theta <- seq(0, pi, length.out = ntheta / 2 + 1)
theta <- c(theta, pi + theta)
theta <- theta[theta != 2*pi]
inc <- 0.00001
theta <- c(theta, inc, pi + inc, pi - inc, 2*pi - inc)
theta <- sort(theta)

coords %<>%
  rowwise() %>%
  do(data_frame(x = .$x, r = .$r, theta = theta)) %>%
  ungroup %>%
  mutate(y = r * cos(theta), z = r * sin(theta))

# get all x & y values used (sort to connect halves on the side)
xs <-
  unique(coords$x) %>%
  sort
ys <-
  unique(coords$y) %>%
  sort

# for each possible x/y pair: get z^2 value
coords <-
  expand.grid(x = xs, y = ys) %>%
  as_data_frame %>%
  mutate(r = r(x), z2 = r^2 - y^2)

# format z coordinates above x/y plane as matrix where columns
# represent x and rows y
zs <- matrix(sqrt(coords$z2), ncol = length(xs), byrow = TRUE)
zs2 <- zs

L <- ncol(zs)
for(i in (L-1):1){
  w <- which(!is.na(zs[, (i+1)]) & is.na(zs[, i]))
  zs2[w, i] <- 0
}

# format x coordiantes as matrix as above (for color gradient)
gradient <-
  rep(xs, length(ys)) %>%
  matrix(ncol = length(xs), byrow = TRUE)

# plot upper half of shape as surface
p <- plot_ly(x = xs, y = ys, z = zs2, surfacecolor = gradient,
             type = "surface", colorbar = list(title = 'x'))

# plot lower have of shape as second surface
p %>%
  add_surface(z = -zs2, showscale = FALSE)

enter image description here

Ответ 4

Это не отвечает на ваш вопрос, но даст результат, с которым вы сможете взаимодействовать на веб-странице: не используйте plot_ly, используйте rgl. Например,

library(rgl)

# Your initial values...

r <- function(x) x^2
int <- c(1, 3)
nx <- 20
ntheta <- 36

# Set up x and colours for each x

x <- seq(int[1], int[2], length.out = nx)
cols <- colorRampPalette(c("blue", "yellow"), space = "Lab")(nx)

clear3d()
shade3d(turn3d(x, r(x), n = ntheta,  smooth = TRUE, 
        material = list(color = rep(cols, each = 4*ntheta))))
aspect3d(1,1,1)
decorate3d()
rglwidget()

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

Вот результат:

enter image description here