使用径向基函数插值球面上的函数

Question

首先,一点背景:

我使用球面谐波作为球体表面上的函数示例,如此图像中的前球体:

我制作了这些球体中的一个,根据其表面上的点处的谐波函数的值着色.我首先做了很多点,所以我的功能非常准确.我称之为我的fine球体.

现在我有了我的fine球体,我在球体上占据了相对较少的点数.这些是我想要插入的点,训练数据,我称之为interp点.以下是我的interp点,在他的值上着色,绘制在我的fine球体上.

现在,该项目的目标是使用这些interp点来训练SciPy径向基函数以在球体上插入我的函数.我能够这样做:

# Train the interpolation using interp coordinates
rbf = Rbf(interp.phi, interp.theta, harmonic13_coarse)
# The result of the interpolation on fine coordinates
interp_values = rbf(fine.phi, fine.theta)

哪个产生了这个插值,绘制在球体上:

希望通过最后一张图片,您可以看到我的问题.注意通过插值的线？这是因为插值数据具有边界.边界是因为我使用球面坐标([0,pi]和[0,2pi]处的边界)训练径向基函数.

rbf = Rbf(interp.phi, interp.theta, harmonic13_coarse)

我的目标,以及为什么我发布这个问题,是使用球体上数据的x,y,z笛卡尔坐标在球体表面上插入我的函数.这样,由于球体没有边界,我不会像在球坐标中那样出现这个边界误差.但是,我只是无法弄清楚如何做到这一点.

我试过简单地给Rbf函数提供x,y,z坐标和函数的值.

rbf=Rbf(interp.x, interp.y, interp.z, harmonic13_coarse)
interp_values=rbf(fine.x,fine.y,fine.z)

但NumPy却给我一个奇异矩阵错误

numpy.linalg.linalg.LinAlgError: singular matrix

有什么方法可以让我用笛卡尔坐标给Rbf我的数据站点,每个站点都有函数值,并且它的行为与球面坐标的行为相似但没有边界吗？从Rbf文档中,有norm用于定义不同距离范数的属性,我是否必须使用球形距离才能使其工作？

我对此很难过.如果您对没有球坐标边界的球体内插函数有任何想法,请告诉我.

这是我的完整代码:

import matplotlib.pyplot as plt
from matplotlib import cm, colors
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from scipy import special
from scipy.interpolate import Rbf
from collections import namedtuple
from mayavi import mlab

# Nice aliases
pi = np.pi
cos = np.cos
sin = np.sin

# Creating a sphere in Cartesian and Sphereical
# Saves coordinates as named tuples


def coordinates(r, n):
    phi, theta = np.mgrid[0:pi:n, 0:2 * pi:n]
    Coor = namedtuple('Coor', 'r phi theta x y z')
    r = r
    x = r * sin(phi) * cos(theta)
    y = r * sin(phi) * sin(theta)
    z = r * cos(phi)
    return Coor(r, phi, theta, x, y, z)

# Creating a sphere
# fine is coordinates on a fine grid
# interp is coordinates on coarse grid for training interpolation
fine = coordinates(1, 100j)
interp = coordinates(1, 5j)


# Defining finection to colour sphere
# Here we are using a spherical harmonic
def harmonic(m, n, theta, phi):
    return special.sph_harm(m, n, theta, phi).real
norm = colors.Normalize()

# One example of the harmonic function, for testing
harmonic13_fine = harmonic(1, 3, fine.theta, fine.phi)
harmonic13_coarse = harmonic(1, 3, interp.theta, interp.phi)


# Train the interpolation using interp coordinates
rbf = Rbf(interp.phi, interp.theta, harmonic13_coarse)
# The result of the interpolation on fine coordinates
interp_values = rbf(fine.phi, fine.theta)


rbf=Rbf(interp.x, interp.y, interp.z, harmonic13_coarse)
interp_values=rbf(fine.x,fine.y,fine.z)

#Figure of harmoinc function on sphere in fine cordinates
#Points3d showing interpolation training points coloured to their value
mlab.figure()
vmax, vmin = np.max(harmonic13_fine), np.min(harmonic13_fine)
mlab.mesh(fine.x, fine.y, fine.z, scalars=harmonic13_fine, vmax=vmax, vmin=vmin)
mlab.points3d(interp.x, interp.y, interp.z, harmonic13_coarse,
              scale_factor=0.1, scale_mode='none', vmax=vmax, vmin=vmin)


#Figure showing results of rbf interpolation
mlab.figure()
vmax, vmin = np.max(harmonic13_fine), np.min(harmonic13_fine)
mlab.mesh(fine.x, fine.y, fine.z, scalars=interp_values)
# mlab.points3d(interp.x, interp.y, interp.z, scalars, scale_factor=0.1, scale_mode='none',vmax=vmax, vmin=vmin)

mlab.show()

Answer 1

您看到的边界是因为您将封闭曲面 (S2) 映射到开放曲面 (R2)。不管怎样，你都会有界限。流形的局部属性是兼容的，因此它适用于大部分球体，但不适用于global，您会得到一条线。

解决方法是使用地图集而不是单个图表。图集是重叠图表的集合。在重叠区域中，您需要定义权重，这是一个在每个图表上从 0 到 1 的平滑函数。（抱歉，可能微分几何不是您期望听到的）。

如果您不想一直走到这里，您会注意到您的原始球体有一个方差最小的赤道。然后，您可以旋转细球体并使其与线重合。它不能解决您的问题，但它肯定可以减轻它。