Python SciPy RectSphereBivariateSpline 插值生成错误数据？

Question

我有一个非常粗糙的球体上的 3D 测量数据，我想对其进行插值。在 stackoverflow 的 @M4rtini 和 @HYRY 的大力帮助下，我现在能够生成工作代码（基于 SciPy 的 RectSphereBivariateSpline 示例中的原始示例）。

测试数据可以在这里找到：testdata

""" read csv input file, post process and plot 3D data """
import csv
import numpy as np
from mayavi import mlab
from scipy.interpolate import RectSphereBivariateSpline

# user input
nElevationPoints = 17 # needs to correspond with csv file
nAzimuthPoints = 40 # needs to correspond with csv file
threshold = - 40 # needs to correspond with how measurement data was captured
turnTableStepSize = 72 # needs to correspond with measurement settings
resolution = 0.125 # needs to correspond with measurement settings

# read data from file
patternData = np.empty([nElevationPoints, nAzimuthPoints]) # empty buffer
ifile  = open('ttest.csv') # need the 'b' suffix to prevent blank rows being inserted
reader = csv.reader(ifile,delimiter=',')
reader.next() # skip first line in csv file as this is only text
for nElevation in range (0,nElevationPoints):
    # azimuth
    for nAzimuth in range(0,nAzimuthPoints):  
        patternData[nElevation,nAzimuth] = reader.next()[2]
ifile.close()

# post process
def r(thetaIndex,phiIndex):
    """r(thetaIndex,phiIndex): function in 3D plotting to return positive vector length from patternData[theta,phi]"""
    radius = -threshold + patternData[thetaIndex,phiIndex]
    return radius

#phi,theta = np.mgrid[0:nAzimuthPoints,0:nElevationPoints]
theta = np.arange(0,nElevationPoints)
phi = np.arange(0,nAzimuthPoints)
thetaMesh, phiMesh = np.meshgrid(theta,phi)
stepSizeRad = turnTableStepSize * resolution * np.pi / 180
theta = theta * stepSizeRad
phi = phi * stepSizeRad

# create new grid to interpolate on
phiIndex = np.arange(1,361)
phiNew = phiIndex*np.pi/180
thetaIndex = np.arange(1,181)
thetaNew = thetaIndex*np.pi/180
thetaNew,phiNew = np.meshgrid(thetaNew,phiNew)
# create interpolator object and interpolate
data = r(thetaMesh,phiMesh)
theta[0] += 1e-6 # zero values for theta cause program to halt; phi makes no sense at theta=0
lut = RectSphereBivariateSpline(theta,phi,data.T)
data_interp = lut.ev(thetaNew.ravel(),phiNew.ravel()).reshape((360,180)).T

def rInterp(theta,phi):
    """rInterp(theta,phi): function in 3D plotting to return positive vector length from interpolated patternData[theta,phi]"""
    thetaIndex = theta/(np.pi/180)
    thetaIndex = thetaIndex.astype(int)
    phiIndex = phi/(np.pi/180)
    phiIndex = phiIndex.astype(int)
    radius = data_interp[thetaIndex,phiIndex]
    return radius
# recreate mesh minus one, needed otherwise the below gives index error, but why??
phiIndex = np.arange(0,360)
phiNew = phiIndex*np.pi/180
thetaIndex = np.arange(0,180)
thetaNew = thetaIndex*np.pi/180
thetaNew,phiNew = np.meshgrid(thetaNew,phiNew)

x = (rInterp(thetaNew,phiNew)*np.cos(phiNew)*np.sin(thetaNew))
y = (-rInterp(thetaNew,phiNew)*np.sin(phiNew)*np.sin(thetaNew))
z = (rInterp(thetaNew,phiNew)*np.cos(thetaNew))

# plot 3D data
obj = mlab.mesh(x, y, z, colormap='jet')
obj.enable_contours = True
obj.contour.filled_contours = True
obj.contour.number_of_contours = 20
mlab.show()

虽然代码运行，但结果图与非插值数据有很大不同，见图

作为参考。

此外，在运行交互式会话时，data_interp 的值 (>3e5) 比原始数据大得多（最大值约为 20）。

有人知道我做错了什么吗？

Answer 1

我好像解决了！

事实上，我试图推断，而我只能对这些分散的数据进行插值。所以新的插值网格应该只上升到 theta = 140 度左右。

但最重要的变化是在 RectSphereBivariateSpline 调用中增加了参数 s=900。

所以我现在有以下代码：

""" read csv input file, post process and plot 3D data """
import csv
import numpy as np
from mayavi import mlab
from scipy.interpolate import RectSphereBivariateSpline

# user input
nElevationPoints = 17 # needs to correspond with csv file
nAzimuthPoints = 40 # needs to correspond with csv file
threshold = - 40 # needs to correspond with how measurement data was captured
turnTableStepSize = 72 # needs to correspond with measurement settings
resolution = 0.125 # needs to correspond with measurement settings

# read data from file
patternData = np.empty([nElevationPoints, nAzimuthPoints]) # empty buffer
ifile  = open('ttest.csv') # need the 'b' suffix to prevent blank rows being inserted
reader = csv.reader(ifile,delimiter=',')
reader.next() # skip first line in csv file as this is only text
for nElevation in range (0,nElevationPoints):
    # azimuth
    for nAzimuth in range(0,nAzimuthPoints):  
        patternData[nElevation,nAzimuth] = reader.next()[2]
ifile.close()

# post process
def r(thetaIndex,phiIndex):
    """r(thetaIndex,phiIndex): function in 3D plotting to return positive vector length from patternData[theta,phi]"""
    radius = -threshold + patternData[thetaIndex,phiIndex]
    return radius

#phi,theta = np.mgrid[0:nAzimuthPoints,0:nElevationPoints]
theta = np.arange(0,nElevationPoints)
phi = np.arange(0,nAzimuthPoints)
thetaMesh, phiMesh = np.meshgrid(theta,phi)
stepSizeRad = turnTableStepSize * resolution * np.pi / 180
theta = theta * stepSizeRad
phi = phi * stepSizeRad

# create new grid to interpolate on
phiIndex = np.arange(1,361)
phiNew = phiIndex*np.pi/180
thetaIndex = np.arange(1,141)
thetaNew = thetaIndex*np.pi/180
thetaNew,phiNew = np.meshgrid(thetaNew,phiNew)
# create interpolator object and interpolate
data = r(thetaMesh,phiMesh)
theta[0] += 1e-6 # zero values for theta cause program to halt; phi makes no sense at theta=0
lut = RectSphereBivariateSpline(theta,phi,data.T,s=900)
data_interp = lut.ev(thetaNew.ravel(),phiNew.ravel()).reshape((360,140)).T

def rInterp(theta,phi):
    """rInterp(theta,phi): function in 3D plotting to return positive vector length from interpolated patternData[theta,phi]"""
    thetaIndex = theta/(np.pi/180)
    thetaIndex = thetaIndex.astype(int)
    phiIndex = phi/(np.pi/180)
    phiIndex = phiIndex.astype(int)
    radius = data_interp[thetaIndex,phiIndex]
    return radius
# recreate mesh minus one, needed otherwise the below gives index error, but why??
phiIndex = np.arange(0,360)
phiNew = phiIndex*np.pi/180
thetaIndex = np.arange(0,140)
thetaNew = thetaIndex*np.pi/180
thetaNew,phiNew = np.meshgrid(thetaNew,phiNew)

x = (rInterp(thetaNew,phiNew)*np.cos(phiNew)*np.sin(thetaNew))
y = (-rInterp(thetaNew,phiNew)*np.sin(phiNew)*np.sin(thetaNew))
z = (rInterp(thetaNew,phiNew)*np.cos(thetaNew))

# plot 3D data
intensity = rInterp(thetaNew,phiNew)
obj = mlab.mesh(x, y, z, scalars = intensity, colormap='jet')
obj.enable_contours = True
obj.contour.filled_contours = True
obj.contour.number_of_contours = 20
mlab.show()

结果图与原始非插值数据很好地比较：

我不完全理解为什么 s 应该设置为 900，因为 RectSphereBivariateSpline 文档说 s=0 用于插值。但是，当进一步阅读文档时，我会有所了解：

选择 s 的最佳值可能是一项微妙的任务。 s 的推荐值取决于数据值的准确性。如果用户对数据的统计误差有所了解，她也可以找到 s 的适当估计值。通过假设，如果她指定正确的 s，则插值器将使用样条 f(u,v) 精确再现数据基础的函数，她可以计算 sum((r(i,j)-s(u(i ),v(j)))**2) 为这个 s 找到一个好的估计。例如，如果她知道她的 r(i,j) 值的统计误差不大于 0.1，她可能会期望一个好的 s 的值不大于 u.size * v.size * (0.1 )**2。 如果对 r(i,j) 中的统计误差一无所知，则 s 必须通过反复试验确定。最好的方法是从一个非常大的 s 值开始（确定最小二乘多项式和 s 的相应上限 fp0），然后逐渐减小 s 的值（例如在开始时减少 10 倍，即s = fp0 / 10, fp0 / 100, ... 更仔细，因为近似值显示更多细节）以获得更紧密的拟合。