如何在球体上绘制等高线？

Question

假设我有一个

rgl

球形物体，我想在该球体的表面上绘制轮廓。我认为

contourLines3d

R 包中的

rgl

函数是我最好的选择，所以我将在这里为我们调整 他们的 MWE：

lat <- matrix(seq(90, -90, len = 50)*pi/180, 50, 50, byrow = TRUE)
long <- matrix(seq(-180, 180, len = 50)*pi/180, 50, 50)

r <- 6378.1 # radius of Earth in km
x <- r*cos(lat)*cos(long)
y <- r*cos(lat)*sin(long)
z <- r*sin(lat)

open3d()
ids <- persp3d(x, y, z, col = "white",
        specular = "black", axes = FALSE, box = FALSE, xlab = "", ylab = "", zlab = "",
        normal_x = x, normal_y = y, normal_z = z, polygon_offset = 1)
        
contourLines3d(ids, func)

我需要正确指定 func 来绘制轮廓。这就是我正在挣扎的地方。假设我有一个 pdf 函数 (

dSphereFunc

) 可以计算球体上任意点的密度。如何使用我的

dSphereFunc

绘制轮廓？我确信有一种相对简单的方法可以做到这一点，但我在这里确实遇到了麻烦。谢谢您的帮助。

更新： 根据评论中的要求，我提供了我正在尝试评估的密度函数：

dSphereFunc = function(y, mu, Vinv){
    d = length(mu)
    ytVinvy = t(y)%*%Vinv%*%y
    ytmu = t(y)%*%mu
    result = ((2*pi)^(-(d-1)/2))/(ytVinvy^(d/2))*exp((ytmu^2/ytVinvy - t(mu)%*%mu)/2)*mdmo(d, ytmu/(ytVinvy^(1/2)))
    return(result)
}
mdmo = function(d, t){
    t = c(t) # formatting from matrix to constant
    integrand = function(x) x^(d-1)*exp(-(x-t)^2/2)
    integral = integrate(integrand, lower = 0, upper = Inf)
    return((2*pi)^(-1/2)*integral$value)
}

在上面的函数中，我们将

mu

（

dx1

向量）和

V

（

dxd

矩阵）视为固定（即已经估计），其中

d=3

;另外，

y

是一个

dx1

向量，它是球坐标 x、y 和 z。

要拥有 MWE，以下是

mu

和

Vinv

的一些值：

mu=c(2,2,2)
Vinv=solve(matrix(c(1.57,-0.08,-0.5,-0.08,0.74,0.34,-0.5,0.34,1.16),nrow=3))
test=c(1,3,4)
test=test/norm(test,type="2")
dSphereFunc(test,mu,Vinv) # 1.156622

我认为我的麻烦是

dSphereFunc

需要三个参数，后两个已修复。因此，我尝试使用以下函数，该函数只需要一个参数，

xyz

，它是球坐标 x、y 和 z 的

dx1

向量：

densFunc = function(xyz){apply(xyz, 1, function(point) dSphereFunc(point, mu, Vinv))}

但是，这仍然不起作用。

Answer 1

您使用的代码很好，唯一的问题是绘图：轮廓全部出现在球体的背面，所以您看不到它们。在绘制球体时设置

alpha = 0.5

可以使它们可用，或者您可以旋转球体直到它们出现。这里的线条从球体的背面显示出来：

lat <- matrix(seq(90, -90, len = 50)*pi/180, 50, 50, byrow = TRUE)
long <- matrix(seq(-180, 180, len = 50)*pi/180, 50, 50)

r <- 6378.1 # radius of Earth in km
x <- r*cos(lat)*cos(long)
y <- r*cos(lat)*sin(long)
z <- r*sin(lat)

library(rgl)
open3d()
#> glX 
#>   1
ids <- persp3d(x, y, z, col = "white",
               specular = "black", axes = FALSE, box = FALSE, xlab = "", ylab = "", zlab = "",
               normal_x = x, normal_y = y, normal_z = z, polygon_offset = 1, alpha = 0.5)

dSphereFunc = function(y, mu, Vinv){
  d = length(mu)
  ytVinvy = t(y)%*%Vinv%*%y
  ytmu = t(y)%*%mu
  result = ((2*pi)^(-(d-1)/2))/(ytVinvy^(d/2))*exp((ytmu^2/ytVinvy - t(mu)%*%mu)/2)*mdmo(d, ytmu/(ytVinvy^(1/2)))
  return(result)
}
mdmo = function(d, t){
  t = c(t) # formatting from matrix to constant
  integrand = function(x) x^(d-1)*exp(-(x-t)^2/2)
  integral = integrate(integrand, lower = 0, upper = Inf)
  return((2*pi)^(-1/2)*integral$value)
}

mu=c(2,2,2)
Vinv=solve(matrix(c(1.57,-0.08,-0.5,-0.08,0.74,0.34,-0.5,0.34,1.16),nrow=3))
test=c(1,3,4)
test=test/norm(test,type="2")
dSphereFunc(test,mu,Vinv) 
#>          [,1]
#> [1,] 1.156622

densFunc = function(xyz){apply(xyz, 1, function(point) dSphereFunc(point, mu, Vinv))}


contourLines3d(ids, densFunc)

^{创建于 2024-05-09，使用 reprex v2.1.0}