如何将等值面夹到球上？

Question

考虑 Togliatti 隐式曲面。我想将它夹在以原点为中心、半径为 4.8 的球上。使用 misc3d 包的解决方案包括使用 computeContour3d 函数的 mask 参数，它只允许使用满足 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))

这是我看到的输出：

它不如你的好，但它适用于许多不同的边界函数。

编辑添加：rgl 的 0.100.26 版现在有一个功能 clipMesh3d，其中包含这些想法。

如何将等值面夹到球上？

How to clip an isosurface to a ball?

r

3d

surface

rgl