用Delaunay三角剖分（n-dim）插值

Question

我想在Python中使用Delaunay三角剖分法來插入3D中的點。

我有什麼是

# my array of points 
points = [[1,2,3], [2,3,4], ...] 
# my array of values 
values = [7, 8, ...] 
# an object with triangulation 
tri = Delaunay(points)   
# a set of points at which I want to interpolate 
p = [[1.5, 2.5, 3.5], ...] 
# this gets simplexes that contain given points 
s = tri.find_simplex(p) 
# this gets vertices for the simplexes 
v = tri.vertices[s] 

我才能夠找到一個answer這裏建議使用transform方法插值，但是沒有被任何更具體。

我需要知道的是如何使用包含單純形的頂點來獲得線性插值的權重。讓我們假設一個一般的n-dim案例，以便答案不依賴於維度。

編輯：我不想使用LinearNDInterpolator或類似的方法，因爲我沒有在每個點的數字作爲一個值，但更復雜的東西（數組/函數）。

Answer 1

一些試驗後，該解決方案看起來很簡單（這post相當有用）：

# dimension of the problem (in this example I use 3D grid, 
# but the method works for any dimension n>=2) 
n = 3 
# my array of grid points (array of n-dimensional coordinates) 
points = [[1,2,3], [2,3,4], ...] 
# each point has some assigned value that will be interpolated 
# (e.g. a float, but it can be a function or anything else) 
values = [7, 8, ...] 
# a set of points at which I want to interpolate (it must be a NumPy array) 
p = np.array([[1.5, 2.5, 3.5], [1.1, 2.2, 3.3], ...]) 

# create an object with triangulation 
tri = Delaunay(points)   
# find simplexes that contain interpolated points 
s = tri.find_simplex(p) 
# get the vertices for each simplex 
v = tri.vertices[s] 
# get transform matrices for each simplex (see explanation bellow) 
m = tri.transform[s] 

# for each interpolated point p, mutliply the transform matrix by 
# vector p-r, where r=m[:,n,:] is one of the simplex vertices to which 
# the matrix m is related to (again, see bellow) 
b = np.einsum('ijk,ik->ij', m[:,:n,:n], p-m[:,n,:]) 

# get the weights for the vertices; `b` contains an n-dimensional vector 
# with weights for all but the last vertices of the simplex 
# (note that for n-D grid, each simplex consists of n+1 vertices); 
# the remaining weight for the last vertex can be copmuted from 
# the condition that sum of weights must be equal to 1 
w = np.c_[b, 1-b.sum(axis=1)] 

的關鍵方法要明白的是transform，這是簡要記錄，但文件說，所有這需要說。對於每個單純形，transform[:,:n,:n]包含變換矩陣，並且transform[:,n,:]包含矩陣與之相關的向量r。看起來，r向量被選爲單純形的最後一個頂點。

另一個棘手的問題是如何讓b，因爲我想要做的是一樣的東西

for i in range(len(p)): b[i] = m[i,:n,:n].dot(p[i]-m[i,n,:]) 

從本質上講，我需要點產品陣列，同時dot給出了兩個陣列的產品。在各個單形環路像上面會的工作，但它可以在一個步驟中完成得更快，對於其中有numpy.einsum：

b = np.einsum('ijk,ik->ij', m[:,:n,:n], p-m[:,n,:]) 

現在，v包含頂點指數各單體和w持有相應的權重。在設定點p獲取插值p_values，我們做的（注：values必須與NumPy陣列此）：

values = np.array(values) 
for i in range(len(p)): p_values[i] = np.inner(values[v[i]], w[i]) 

或者，我們可以使用`np.einsum」再這樣做在一個單一的步驟：

p_values = np.einsum('ij,ij->i', values[v], w) 

當某些插值點位於網格之外時，必須注意一些情況。在這種情況下，find_simplex(p)返回-1這些點，然後你將不得不掩蓋它們（也許使用masked arrays）。

Answer 2

您需要一個區間和一個線性插值，即邊的長度和插值點距開始頂點的距離。

Answer 3

你並不需要從頭實現這一點，已經有內置的支持SciPy的這個特性：

scipy.interpolate.LinearNDInterpolator