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

如何解决在 3D 球体上插入非均匀分布的点

我在单位球面上有几个点，它们根据 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,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.0001 到 100000，上面的代码总是引发上述异常或：

  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,r_interp)

# Creating the finer theta and phi values for the final plot
plot_points = 100
theta = numpy.linspace(0,plot_points,numpy.ones(plot_points))

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()

这段代码运行流畅，结果如下：

插值看起来没问题除了在不连续的一行上，就像在引导我上这门课的问题中一样。 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,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,z,r_interp)

# predict
heatmap = interpolator(X,Y,Z)

结果如下：

ax.plot_surface(
    X,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)

为了更好地看到差异， 我堆叠了两个图像，减去它们并反转图层。