首先,介绍一下背景:
我使用球谐函数作为球体表面上的函数示例,如下图中的前球体所示:
![enter image description here](https://i.stack.imgur.com/gJLxn.jpg)
我制作了其中一个球体,根据其表面各点的调和函数值进行着色。我首先对大量点执行此操作,因此我的函数非常准确。我把这称为我的fine
领域。
![enter image description here](https://i.stack.imgur.com/fUWbU.png)
现在我有了我的fine
球体,我在球体上取相对少量的点。这些是我希望从训练数据中插入的点,我将它们称为interp
点。这是我的interp
点,根据它们的值着色,绘制在我的fine
sphere.
![enter image description here](https://i.stack.imgur.com/pMHh5.png)
现在,该项目的目标是使用这些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)
Which produced this interpolation, plotted on the sphere:
![enter image description here](https://i.stack.imgur.com/Y92SV.png)
希望通过最后一张图片,您可以看到我的问题。注意到穿过插值的线了吗?这是因为插值数据有边界。边界是因为我使用球坐标(边界在 [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()