如何在Python中为某些球面坐标的点生成高斯分布的偏移新点

Question

我正在处理一些球坐标系中的点。我需要生成新的点作为它们的误差点，以及对旧点的一种偏移。新点应该与旧点保持特定的距离，这个距离是根据高斯分布来分配的。新点与旧点之间的角度并不重要。我想为r方向生成新点，不管phi和theta是什么（球坐标系）。

为了生成根据高斯函数分布的新点，我尝试使用numpy.rand.normal(mean,std,..)。但是它生成的是在均值附近的1D随机点，而这个均值是一个实数。在我的情况下，我需要一种方法来指定旧点的位置，并且我有一个给定的标准差来表示与原始点的距离。

老实说，我没有我代码的副本。它在大学的服务器上。但我们假设我有一个大小为100*3的数组，其中包含一些在圆柱体表面上的点的球坐标（或笛卡尔坐标）。在球坐标的情况下，第一列表示半径值，第二列是theta，第三列显示点的phi。现在我想使用高斯分布从这些点中生成随机点。高斯分布有一个给定的标准差。唯一重要的是，生成的新点在r值上是有限制的。无论点在theta和phi方面的位置如何。

当我尝试numpy.rand.normal(mean,std,..)时，它生成了一些在均值附近的随机点。这对我没有帮助。我想要的是在给定标准差的情况下，生成在旧点上的新点。

任何想法都将不胜感激。

这是一个与我类似的代码，由Ophion编写，链接在这里：如何在圆柱面上生成规则点

def make_cylinder(radius, length, nlength, alpha, nalpha, center, orientation):

#Create the length array
I = np.linspace(0, length, nlength)

#Create alpha array avoid duplication of endpoints
#Conditional should be changed to meet your requirements
if int(alpha) == 360:
    A = np.linspace(0, alpha, num=nalpha, endpoint=False)/180*np.pi
else:
    A = np.linspace(0, alpha, num=nalpha)/180*np.pi

#Calculate X and Y
X = radius * np.cos(A)
Y = radius * np.sin(A)

#Tile/repeat indices so all unique pairs are present
pz = np.tile(I, nalpha)
px = np.repeat(X, nlength)
py = np.repeat(Y, nlength)

points = np.vstack(( pz, px, py )).T

#Shift to center
shift = np.array(center) - np.mean(points, axis=0)
points += shift

#Orient tube to new vector

#Grabbed from an old unutbu answer
def rotation_matrix(axis,theta):
    a = np.cos(theta/2)
    b,c,d = -axis*np.sin(theta/2)
    return np.array([[a*a+b*b-c*c-d*d, 2*(b*c-a*d), 2*(b*d+a*c)],
                     [2*(b*c+a*d), a*a+c*c-b*b-d*d, 2*(c*d-a*b)],
                     [2*(b*d-a*c), 2*(c*d+a*b), a*a+d*d-b*b-c*c]])

ovec = orientation / np.linalg.norm(orientation)
cylvec = np.array([1,0,0])

if np.allclose(cylvec, ovec):
    return points

#Get orthogonal axis and rotation
oaxis = np.cross(ovec, cylvec)
rot = np.arccos(np.dot(ovec, cylvec))

R = rotation_matrix(oaxis, rot)
return points.dot(R)

现在调用这个函数：

points = make_cylinder(3, 5, 5, 360, 10, [0,2,0], [1,0,0])

sigma = 0.5 # 给定的标准差 ossfet_points = numpy.random.normal(np.mean(point[:,0]), sigma, size=(n,3))

Answer 1

如果我没理解错的话，你想要在一个球面上生成随机点，并且这些点距离中心的分布是高斯分布。要做到这一点，你可以使用 numpy.rand.normal 来生成半径的高斯值，这样就解决了后面的问题。

生成随机的球面点稍微复杂一点，不过这里有一些代码可以帮你实现（还有一些数学背景的描述可以在 Wolfram MathWorld 找到）：

import numpy as np

num_points = 500
U = np.random.random(num_points)
V = np.random.random(num_points)

import math as m

def spherical_to_cartesian(vec):                                                  
    '''
    Convert spherical polar coordinates to cartesian coordinates:                       

    See the definition of spherical_cartesian_to_polar.                                 

    @param vec: A vector of the 3 polar coordinates (r, u, v)                           
    @return: (x, y, z)                                                                  
    ''' 
    (r, u, v) = vec                                                                     

    x = r * m.sin(u) * m.cos(v)                                                         
    y = r * m.sin(u) * m.sin(v)                                                         
    z = r * m.cos(u)                                                                    

    return [x, y, z]  

radius = 1.
points = np.array([spherical_to_cartesian([radius, 2 * np.pi * u, np.arccos(2*v - 1)]) for u,v in zip(U,V)])

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax = Axes3D(fig)

ax.plot(points[:,0], points[:,1], points[:,2], 'o')

这样你就能得到像这样的点：

如果你想让这些点的半径符合正态分布，只需要在使用半径变量的列表推导中替换你随机生成的值，像这样：

radii = np.random.normal(10, 3, 100)
points = np.array([spherical_to_cartesian([r, 2 * np.pi * u, np.arccos(2*v - 1)]) for r,u,v in zip(radii, U,V)])

这大概就是你想要的效果吗？