Python SciPy RectSphereBivariateSpline插值生成錯誤的數據？

Question

我有一個非常粗糙的球體上的3D測量數據，我想插入。有了@ M4rtini和@HYRY在stackoverflow的幫助，我現在已經能夠生成工作代碼（基於來自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 max）大得多。

有沒有人知道我可能做錯了什麼？

Answer 1

我似乎已經解決了它！

對於事物，我試圖外推，而我只能插入這些分散的數據。所以新的插值網格應該只能達到θ= 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 ），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，...以及更仔細的近似顯示更多細節）以獲得更接近的擬合。