在R中的特定球体坐标内随机生成点

Question

我有特定的 x,y,z 坐标。我想在给定 x 作为中心和从另一个数据帧 x2 作为半径边缘的球体中生成随机点（因此从 x 到 x2 的距离将是球体半径的长度）。

我已经看到很多关于如何在数学上适当地做到这一点的讨论（随机分布点以避免聚集），并且能够编译最简单的示例 here 和 here 以获取示例 R 代码。

我还发现了这个 [R package sphereplot] (https://cran.r-project.org/web/packages/sphereplot/sphereplot.pdf)，它可能更容易，但很难理解如何应用它。

这些都是很好的起点，但是使用下面的示例代码我不确定如何将其应用于特定的起点/球坐标？

set.seed(101)
n <- 50
theta <- runif(n,0,2*pi)
u <- runif(n,-1,1)
x <- sqrt(1-u^2)*cos(theta)
y <- sqrt(1-u^2)*sin(theta)
z <- u

仅使用我的数据框中的一组/行 x、y、z 坐标：

x = -0.0684486861
y= 0.0125857380
z= 0.0201056441

x2= -0.0684486861
y2 = 0.0125857380
z2= -0.0228805516

我希望 x,y,z 是球体的中心，到 x2,y2,z2 的距离是球体的半径长度/边缘。然后从球体内以 x,y,z 为中心生成随机点。

最终，我尝试使用 100 个球体来比较第二组坐标中的所有点是否在空间中以相似的角度/方向移动。

感谢您的指导。

Answer 1

好吧，让我们将问题拆分为几个子问题。

首先，生成均匀分布在球体上（体积上或表面上）的点，中心位于 (0,0,0) 并具有给定的半径。跟随http://mathworld.wolfram.com/SpherePointPicking.html，与您显示的代码非常接近，

rsphere <- function(n, r = 1.0, surface_only = FALSE) {
    phi       <- runif(n, 0.0, 2.0 * pi)
    cos_theta <- runif(n, -1.0, 1.0)
    sin_theta <- sqrt((1.0-cos_theta)*(1.0+cos_theta))
    radius <- r
    if (surface_only == FALSE) {
        radius <- r * runif(n, 0.0, 1.0)^(1.0/3.0)
    }

    x <- radius * sin_theta * cos(phi)
    y <- radius * sin_theta * sin(phi)
    z <- radius * cos_theta

    cbind(x, y, z)
}

set.seed(312345)
sphere_points <- rsphere(10000)

第二个问题 - 将这些点移动到点 X 的中心

rsphere <- function(n, r = 1.0, surface_only = FALSE, center=cbind(Xx, Xy, Xz)) {
    ....
    cbind(x+center[1], y+center[2], z+center[3])
}

第三个问题 - 在 (Xx, Xy, Xz) 和曲面点 (Yx, Yy, Yz) 处计算半径给定的中心）

radius <- sqrt((Xx-Yx)**2+(Xy-Yy)**2+(Xz-Yz)**2)

将它们组合在一起以获得完全的满足感。好的，现在您提供了中心和半径的值，让我们把它们放在一起

rsphere <- function(n, r = 1.0, surface_only = FALSE, center=cbind(0.0, 0.0, 0.0)) {
    phi       <- runif(n, 0.0, 2.0 * pi)
    cos_theta <- runif(n, -1.0, 1.0)
    sin_theta <- sqrt((1.0-cos_theta)*(1.0+cos_theta))
    radius <- r
    if (surface_only == FALSE) {
        radius <- r * runif(n, 0.0, 1.0)^(1.0/3.0)
    }

    x <- radius * sin_theta * cos(phi)
    y <- radius * sin_theta * sin(phi)
    z <- radius * cos_theta

    # if radius is fixed, we could check it
    # rr = sqrt(x^2+y^2+z^2)
    # print(rr)

    cbind(x+center[1], y+center[2], z+center[3])
}

x1 = -0.0684486861
y1 = 0.0125857380
z1 = 0.0201056441

x2 = -0.0684486861
y2 = 0.0125857380
z2 = -0.0228805516

R = sqrt((x2-x1)^2 + (y2-y1)^2 + (z2-z1)^2)
print(R)

set.seed(32345)
sphere_points <- rsphere(100000, R, FALSE, cbind(x1, y1, z1))

它看起来怎么样？

更新

在表面和体积中各生成 10 个点并打印出来，radius=2 对我来说看起来不错

# 10 points uniform on surface, supposed to have fixed radius
sphere_points <- rsphere(10, 2, TRUE, cbind(x1, y1, z1))
for (k in 1:10) {
    rr <- sqrt((sphere_points[k,1]-x1)^2+(sphere_points[k,2]-y1)^2+(sphere_points[k,3]-z1)^2)
    print(rr)
}

# 10 points uniform in the sphere, supposed to have varying radius
sphere_points <- rsphere(10, 2, FALSE, cbind(x1, y1, z1))
for (k in 1:10) {
    rr <- sqrt((sphere_points[k,1]-x1)^2+(sphere_points[k,2]-y1)^2+(sphere_points[k,3]-z1)^2)
    print(rr)
}

得到

[1] 2
[1] 2
[1] 2
[1] 2
[1] 2
[1] 2
[1] 2
[1] 2
[1] 2
[1] 2

和

[1] 1.32571
[1] 1.505066
[1] 1.255023
[1] 1.82773
[1] 1.219957
[1] 1.641258
[1] 1.881937
[1] 1.083975
[1] 0.4745712
[1] 1.900066

你得到了什么？