Как прикрепить изоповерхность к мячу?

Question

Рассмотрим неявную поверхность Тольятти. Я хочу закрепить его на шаре с центром в начале координат радиусом 4.8. Решение с пакетом misc3d состоит в использовании аргумента mask функции computeContour3d, который позволяет использовать только точки, удовлетворяющие x^2+y^2+z^2 < 4.8^2:

library(misc3d)

# Togliatti surface equation: f(x,y,z) = 0
f <- function(x,y,z){
  w <- 1
  64*(x-w)*
    (x^4-4*x^3*w-10*x^2*y^2-4*x^2*w^2+16*x*w^3-20*x*y^2*w+5*y^4+16*w^4-20*y^2*w^2) - 
    5*sqrt(5-sqrt(5))*(2*z-sqrt(5-sqrt(5))*w)*(4*(x^2+y^2-z^2)+(1+3*sqrt(5))*w^2)^2
}

# make grid
nx <- 220; ny <- 220; nz <- 220
x <- seq(-5, 5, length=nx) 
y <- seq(-5, 5, length=ny)
z <- seq(-4, 4, length=nz) 
g <- expand.grid(x=x, y=y, z=z)

# calculate voxel
voxel <- array(with(g, f(x,y,z)), dim = c(nx,ny,nz))

# mask: keep points satisfying x^2+y^2+z^2 < 4.8^2, in order to 
#       clip the surface to the ball of radius 4.8
mask <- array(with(g, x^2+y^2+z^2 < 4.8^2), dim = c(nx,ny,nz))

# compute isosurface
surf <- computeContour3d(voxel, maxvol=max(voxel), level=0, mask=mask, x=x, y=y, z=z)

# draw isosurface
drawScene.rgl(makeTriangles(surf, smooth=TRUE))

Но границы получаемой поверхности нерегулярны:

Как получить правильные, гладкие границы?

Answer 1

Решение, которое я нашел, обращается к сферическим координатам . Он состоит в определении функции f в терминах сферических координат (ρ, θ, ϕ), затем для вычисления изоповерхности с ρ от 0 до желаемого радиуса, а затем для преобразования результата в декартовы координаты:

# Togliatti surface equation with spherical coordinates
f <- function(ρ, θ, ϕ){
  w <- 1
  x <- ρ*cos(θ)*sin(ϕ)
  y <- ρ*sin(θ)*sin(ϕ)
  z <- ρ*cos(ϕ)
  64*(x-w)*
    (x^4-4*x^3*w-10*x^2*y^2-4*x^2*w^2+16*x*w^3-20*x*y^2*w+5*y^4+16*w^4-20*y^2*w^2) - 
    5*sqrt(5-sqrt(5))*(2*z-sqrt(5-sqrt(5))*w)*(4*(x^2+y^2-z^2)+(1+3*sqrt(5))*w^2)^2
}

# make grid
nρ <- 300; nθ <- 400; nϕ <- 300
ρ <- seq(0, 4.8, length = nρ) # ρ runs from 0 to the desired radius
θ <- seq(0, 2*pi, length = nθ)
ϕ <- seq(0, pi, length = nϕ) 
g <- expand.grid(ρ=ρ, θ=θ, ϕ=ϕ)

# calculate voxel
voxel <- array(with(g, f(ρ,θ,ϕ)), dim = c(nρ,nθ,nϕ))

# calculate isosurface
surf <- computeContour3d(voxel, maxvol=max(voxel), level=0, x=ρ, y=θ, z=ϕ)

# transform to Cartesian coordinates
surf <- t(apply(surf, 1, function(rtp){
  ρ <- rtp[1]; θ <- rtp[2]; ϕ <- rtp[3] 
  c(
    ρ*cos(θ)*sin(ϕ),
    ρ*sin(θ)*sin(ϕ),
    ρ*cos(ϕ)
  )
}))

# draw isosurface
drawScene.rgl(makeTriangles(surf, smooth=TRUE, color = "violetred"))

Теперь полученная поверхность имеет ровные, гладкие границы:

Answer 2

Ваше решение отлично подходит для поставленной вами задачи, потому что сферические координаты настолько естественны для этой границы.Однако здесь есть более общее решение, которое будет работать для других гладких границ.

Идея состоит в том, чтобы разрешить ввод граничной функции и точек отбора, когда они слишком велики или слишком малы.В вашем случае это будет квадратное расстояние от начала координат, и вы захотите отбирать точки, где значение больше 4,8 ^ 2.Но иногда треугольники, нарисованные для получения гладкой поверхности, следует отбирать только частично: одна точка будет сохранена, а две удалены, или две сохранены и одна удалена.Если вы отбросите весь треугольник, который ведет к неровным краям на исходном графике.

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

Эта функция делает это, предполагая, что ввод surf является выводом computeContour3d:

boundSurface <- function(surf, boundFn, bound = 0, greater = TRUE) {
  # Surf is n x 3:  each row is a point, triplets are triangles
  values <- matrix(boundFn(surf) - bound, 3)
  # values is (m = n/3) x 3:  each row is the boundFn value at one point
  # of a triangle
  if (!greater) 
    values <- -values
  keep <- values >= 0
  # counts is m vector counting number of points to keep in each triangle
  counts <- apply(keep, 2, sum)
  # result is initialized to an empty array
  result <- matrix(nrow = 0, ncol = 3)
  # singles is set to all the rows of surf where exactly one
  # point in the triangle is kept, say s x 3
  singles <- surf[rep(counts == 1, each = 3),]
  if (length(singles)) {
    # singleValues is a subset of values where only one vertex is kept
    singleValues <- values[, counts == 1]
    singleIndex <- 3*col(singleValues) + 1:3 - 3
    # good is the index of the vertex to keep, bad are those to fix
    good <- apply(singleValues, 2, function(col) which(col >= 0))
    bad <- apply(singleValues, 2, function(col) which(col < 0))
    for (j in 1:ncol(singleValues)) {
      goodval <- singleValues[good[j], j]
      for (i in 1:2) {
        badval <- singleValues[bad[i,j], j]
        alpha <- goodval/(goodval - badval)
        singles[singleIndex[bad[i,j], j], ] <- 
          (1-alpha)*singles[singleIndex[good[j], j],] +
             alpha *singles[singleIndex[bad[i,j], j],]
      }
    }
    result <- rbind(result, singles)
  }
  doubles <- surf[rep(counts == 2, each = 3),]
  if (length(doubles)) {
    # doubleValues is a subset of values where two vertices are kept
    doubleValues <- values[, counts == 2]
    doubleIndex <- 3*col(doubleValues) + 1:3 - 3
    doubles2 <- doubles
    # good is the index of the vertex to keep, bad are those to fix
    good <- apply(doubleValues, 2, function(col) which(col >= 0))
    bad <- apply(doubleValues, 2, function(col) which(col < 0))
    newvert <- matrix(NA, 2, 3)
    for (j in 1:ncol(doubleValues)) {
      badval <- doubleValues[bad[j], j]
      for (i in 1:2) {
        goodval <- doubleValues[good[i,j], j]
        alpha <- goodval/(goodval - badval)
        newvert[i,] <- 
          (1-alpha)*doubles[doubleIndex[good[i,j], j],] +
             alpha *doubles[doubleIndex[bad[j], j],]
      }
      doubles[doubleIndex[bad[j], j],] <- newvert[1,]
      doubles2[doubleIndex[good[1,j], j],] <- newvert[1,]
      doubles2[doubleIndex[bad[j], j],] <- newvert[2,]
    }
    result <- rbind(result, doubles, doubles2)
  }
  # Finally add all the rows of surf where the whole
  # triangle is kept
  rbind(result, surf[rep(counts == 3, each = 3),])
}

Вы будете использовать его после computeContour3d и до makeTriangles, например,

fn <- function(x) { 
  apply(x^2, 1, sum)
}

drawScene.rgl(makeTriangles(boundSurface(surf, fn, bound = 4.8^2, 
                                         greater = FALSE), 
                            smooth = TRUE))

Вот вывод, который я вижу:

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

Отредактировано для добавления : версия 0.100.26 изrgl теперь имеет функцию clipMesh3d, которая включает эти идеи.