gpt4 book ai didi

python - 在 3D 球体上插入非均匀分布的点

转载 作者:行者123 更新时间:2023-12-04 11:35:05 26 4
gpt4 key购买 nike

我在单位球面上有几个点根据 https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf 中描述的算法分布(并在下面的代码中实现)。在这些点中的每一个上,我都有一个值,在我的特定情况下,它表示 1 减去一个小错误。错误在 [0, 0.1]如果这很重要,那么我的值在 [0.9, 1] .
可悲的是,计算错误是一个代价高昂的过程,我无法根据需要计算尽可能多的点。不过,我希望我的情节看起来像我在绘制“连续”的东西。
所以我想为我的数据拟合一个插值函数,以便能够根据需要采样尽可能多的点。
经过一番研究,我发现 scipy.interpolate.SmoothSphereBivariateSpline这似乎正是我想要的。但我不能让它正常工作。
问题:我可以用什么来插值(样条,线性插值,目前什么都可以)我在单位球体上的数据?答案可以是“您误用了 scipy.interpolation ,这是执行此操作的正确方法”或“此其他功能更适合您的问题”。
应该可以使用 numpy 执行的示例代码和 scipy安装:

import typing as ty

import numpy
import scipy.interpolate


def get_equidistant_points(N: int) -> ty.List[numpy.ndarray]:
"""Generate approximately n points evenly distributed accros the 3-d sphere.

This function tries to find approximately n points (might be a little less
or more) that are evenly distributed accros the 3-dimensional unit sphere.

The algorithm used is described in
https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf.
"""
# Unit sphere
r = 1

points: ty.List[numpy.ndarray] = list()

a = 4 * numpy.pi * r ** 2 / N
d = numpy.sqrt(a)
m_v = int(numpy.round(numpy.pi / d))
d_v = numpy.pi / m_v
d_phi = a / d_v

for m in range(m_v):
v = numpy.pi * (m + 0.5) / m_v
m_phi = int(numpy.round(2 * numpy.pi * numpy.sin(v) / d_phi))
for n in range(m_phi):
phi = 2 * numpy.pi * n / m_phi
points.append(
numpy.array(
[
numpy.sin(v) * numpy.cos(phi),
numpy.sin(v) * numpy.sin(phi),
numpy.cos(v),
]
)
)
return points


def cartesian2spherical(x: float, y: float, z: float) -> numpy.ndarray:
r = numpy.linalg.norm([x, y, z])
theta = numpy.arccos(z / r)
phi = numpy.arctan2(y, x)
return numpy.array([r, theta, phi])


n = 100
points = get_equidistant_points(n)
# Random here, but costly in real life.
errors = numpy.random.rand(len(points)) / 10

# Change everything to spherical to use the interpolator from scipy.
ideal_spherical_points = numpy.array([cartesian2spherical(*point) for point in points])
r_interp = 1 - errors
theta_interp = ideal_spherical_points[:, 1]
phi_interp = ideal_spherical_points[:, 2]
# Change phi coordinate from [-pi, pi] to [0, 2pi] to please scipy.
phi_interp[phi_interp < 0] += 2 * numpy.pi

# Create the interpolator.
interpolator = scipy.interpolate.SmoothSphereBivariateSpline(
theta_interp, phi_interp, r_interp
)

# Creating the finer theta and phi values for the final plot
theta = numpy.linspace(0, numpy.pi, 100, endpoint=True)
phi = numpy.linspace(0, numpy.pi * 2, 100, endpoint=True)

# Creating the coordinate grid for the unit sphere.
X = numpy.outer(numpy.sin(theta), numpy.cos(phi))
Y = numpy.outer(numpy.sin(theta), numpy.sin(phi))
Z = numpy.outer(numpy.cos(theta), numpy.ones(100))

thetas, phis = numpy.meshgrid(theta, phi)
heatmap = interpolator(thetas, phis)
上面代码的问题:
  • 使用原样的代码,我有一个
    ValueError: The required storage space exceeds the available storage space: nxest or nyest too small, or s too small. The weighted least-squares spline corresponds to the current set of knots.
    在初始化 interpolator 时引发实例。
  • 上面的问题好像说我应该改变s的值那是 scipy.interpolate.SmoothSphereBivariateSpline 的参数之一.我测试了 s 的不同值范围从 0.0001100000 , 上面的代码总是引发上述异常或:
    ValueError: Error code returned by bispev: 10

  • 编辑:我在这里包括我的发现。它们不能真正被视为解决方案,这就是为什么我正在编辑而不是作为答案发布。
    通过更多的研究,我发现了这个问题 Using Radial Basis Functions to Interpolate a Function on a Sphere .作者和我有完全一样的问题,使用了不同的插值器: scipy.interpolate.Rbf .我通过替换插值器和绘图更改了上面的代码:
    # Create the interpolator.
    interpolator = scipy.interpolate.Rbf(theta_interp, phi_interp, r_interp)

    # Creating the finer theta and phi values for the final plot
    plot_points = 100
    theta = numpy.linspace(0, numpy.pi, plot_points, endpoint=True)
    phi = numpy.linspace(0, numpy.pi * 2, plot_points, endpoint=True)

    # Creating the coordinate grid for the unit sphere.
    X = numpy.outer(numpy.sin(theta), numpy.cos(phi))
    Y = numpy.outer(numpy.sin(theta), numpy.sin(phi))
    Z = numpy.outer(numpy.cos(theta), numpy.ones(plot_points))

    thetas, phis = numpy.meshgrid(theta, phi)
    heatmap = interpolator(thetas, phis)


    import matplotlib as mpl
    import matplotlib.pyplot as plt
    from matplotlib import cm

    colormap = cm.inferno
    normaliser = mpl.colors.Normalize(vmin=numpy.min(heatmap), vmax=1)
    scalar_mappable = cm.ScalarMappable(cmap=colormap, norm=normaliser)
    scalar_mappable.set_array([])

    fig = plt.figure()
    ax = fig.add_subplot(111, projection="3d")
    ax.plot_surface(
    X,
    Y,
    Z,
    facecolors=colormap(normaliser(heatmap)),
    alpha=0.7,
    cmap=colormap,
    )
    plt.colorbar(scalar_mappable)
    plt.show()
    此代码运行顺利并给出以下结果:
    enter image description here
    插值似乎没问题 除了 在一条不连续的线上,就像在引导我上这门课的问题中一样。 One of the answer给出使用不同距离的想法,更适应球面坐标:Haversine 距离。
    def haversine(x1, x2):
    theta1, phi1 = x1
    theta2, phi2 = x2
    return 2 * numpy.arcsin(
    numpy.sqrt(
    numpy.sin((theta2 - theta1) / 2) ** 2
    + numpy.cos(theta1) * numpy.cos(theta2) * numpy.sin((phi2 - phi1) / 2) ** 2
    )
    )


    # Create the interpolator.
    interpolator = scipy.interpolate.Rbf(theta_interp, phi_interp, r_interp, norm=haversine)
    执行时会发出警告:
    LinAlgWarning: Ill-conditioned matrix (rcond=1.33262e-19): result may not be accurate.
    self.nodes = linalg.solve(self.A, self.di)
    并且结果完全不是预期的:插值函数的值可能高达 -1这显然是错误的。

    最佳答案

    您可以使用 笛卡尔坐标而不是球坐标。
    Rbf 使用的默认范数参数 ( 'euclidean' )足够了

    # interpolation
    x, y, z = numpy.array(points).T
    interpolator = scipy.interpolate.Rbf(x, y, z, r_interp)

    # predict
    heatmap = interpolator(X, Y, Z)
    结果如下:
    ax.plot_surface(
    X, Y, Z,
    rstride=1, cstride=1,
    # or rcount=50, ccount=50,
    facecolors=colormap(normaliser(heatmap)),
    cmap=colormap,
    alpha=0.7, shade=False
    )
    ax.set_xlabel('x axis')
    ax.set_ylabel('y axis')
    ax.set_zlabel('z axis')


    如果需要,您还可以使用余弦距离(范数参数):
    def cosine(XA, XB):
    if XA.ndim == 1:
    XA = numpy.expand_dims(XA, axis=0)
    if XB.ndim == 1:
    XB = numpy.expand_dims(XB, axis=0)
    return scipy.spatial.distance.cosine(XA, XB)
    Cosine Distance
    为了更好地看到差异,
    我堆叠了两个图像,减去它们并反转图层。
    enter image description here

    关于python - 在 3D 球体上插入非均匀分布的点,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/68199445/

    26 4 0
    Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
    广告合作:1813099741@qq.com 6ren.com