如果我想要随机非结构化数据的 3D 样条/平滑插值怎么办？

Posted 2023-03-12

【中文标题】如果我想要随机非结构化数据的 3D 样条/平滑插值怎么办？

【英文标题】：What to do if I want 3D spline/smooth interpolation of random unstructured data?

【发布时间】：2015-12-21 14:09:58

【问题描述】：

我受到@James 的answer 的启发，想看看如何使用griddata 和map_coordinates。在下面的示例中，我展示的是 2D 数据，但我的兴趣是 3D。我注意到griddata 仅提供 1D 和 2D 的样条曲线，并且仅限于 3D 及更高版本的线性插值（可能有很好的理由）。但是，map_coordinates 似乎适用于使用高阶（比分段线性更平滑）插值的 3D。

我的主要问题：如果我在 3D 中有随机的非结构化数据（我不能使用 map_coordinates），有没有办法在 NumPy SciPy 宇宙中比分段线性插值更平滑，或者至少在附近？

我的第二个问题：griddata 中没有用于 3D 的样条曲线是因为实施起来困难或乏味，还是存在根本性的困难？

下面的图像和可怕的 python 显示了我目前对 griddata 和 map_coordinates 如何使用或不能使用的理解。沿着粗黑线进行插值。

结构化数据：

联合国结构化数据：

可怕的蟒蛇：

import numpy as np
import matplotlib.pyplot as plt

def g(x, y):
    return np.exp(-((x-1.0)**2 + (y-1.0)**2))

def findit(x, X):  # or could use some 1D interpolation
    fraction = (x - X[0]) / (X[-1]-X[0])
    return fraction * float(X.shape[0]-1)

nth, nr = 12, 11
theta_min, theta_max = 0.2, 1.3
r_min,     r_max     = 0.7, 2.0

theta = np.linspace(theta_min, theta_max, nth)
r     = np.linspace(r_min,     r_max,     nr)

R, TH = np.meshgrid(r, theta)
Xp, Yp  = R*np.cos(TH), R*np.sin(TH)
array = g(Xp, Yp)

x, y = np.linspace(0.0, 2.0, 200), np.linspace(0.0, 2.0, 200)
X, Y = np.meshgrid(x, y)
blob = g(X, Y)

xtest = np.linspace(0.25, 1.75, 40)
ytest = np.zeros_like(xtest) + 0.75

rtest = np.sqrt(xtest**2 + ytest**2)
thetatest = np.arctan2(xtest, ytest)

ir = findit(rtest, r)
it = findit(thetatest, theta)

plt.figure()

plt.subplot(2,1,1)

plt.scatter(100.0*Xp.flatten(), 100.0*Yp.flatten())
plt.plot(100.0*xtest, 100.0*ytest, '-k', linewidth=3)
plt.hold
plt.imshow(blob, origin='lower', cmap='gray')

plt.text(5, 5, "don't use jet!", color='white')


exact = g(xtest, ytest)

import scipy.ndimage.interpolation as spndint
ndint0 = spndint.map_coordinates(array, [it, ir], order=0)
ndint1 = spndint.map_coordinates(array, [it, ir], order=1)
ndint2 = spndint.map_coordinates(array, [it, ir], order=2)

import scipy.interpolate as spint
points = np.vstack((Xp.flatten(), Yp.flatten())).T   # could use np.array(zip(...))
grid_x = xtest
grid_y = np.array([0.75])

g0 = spint.griddata(points, array.flatten(), (grid_x, grid_y), method='nearest')
g1 = spint.griddata(points, array.flatten(), (grid_x, grid_y), method='linear')
g2 = spint.griddata(points, array.flatten(), (grid_x, grid_y), method='cubic')


plt.subplot(4,2,5)

plt.plot(exact, 'or')
#plt.plot(ndint0)
plt.plot(ndint1)
plt.plot(ndint2)
plt.title("map_coordinates")

plt.subplot(4,2,6)

plt.plot(exact, 'or')
#plt.plot(g0)
plt.plot(g1)
plt.plot(g2)
plt.title("griddata")

plt.subplot(4,2,7)

#plt.plot(ndint0 - exact)
plt.plot(ndint1 - exact)
plt.plot(ndint2 - exact)
plt.title("error map_coordinates")

plt.subplot(4,2,8)

#plt.plot(g0 - exact)
plt.plot(g1 - exact)
plt.plot(g2 - exact)
plt.title("error griddata")

plt.show()


seed_points_rand = 2.0 * np.random.random((400, 2))
rr = np.sqrt((seed_points_rand**2).sum(axis=-1))
thth = np.arctan2(seed_points_rand[...,1], seed_points_rand[...,0])
isinside = (rr>r_min) * (rr<r_max) * (thth>theta_min) * (thth<theta_max)
points_rand = seed_points_rand[isinside]

Xprand, Yprand = points_rand.T  # unpack
array_rand = g(Xprand, Yprand)

grid_x = xtest
grid_y = np.array([0.75])

plt.figure()

plt.subplot(2,1,1)

plt.scatter(100.0*Xprand.flatten(), 100.0*Yprand.flatten())
plt.plot(100.0*xtest, 100.0*ytest, '-k', linewidth=3)
plt.hold
plt.imshow(blob, origin='lower', cmap='gray')
plt.text(5, 5, "don't use jet!", color='white')


g0rand = spint.griddata(points_rand, array_rand.flatten(), (grid_x, grid_y), method='nearest')
g1rand = spint.griddata(points_rand, array_rand.flatten(), (grid_x, grid_y), method='linear')
g2rand = spint.griddata(points_rand, array_rand.flatten(), (grid_x, grid_y), method='cubic')

plt.subplot(4,2,6)

plt.plot(exact, 'or')
#plt.plot(g0rand)
plt.plot(g1rand)
plt.plot(g2rand)
plt.title("griddata")


plt.subplot(4,2,8)

#plt.plot(g0rand - exact)
plt.plot(g1rand - exact)
plt.plot(g2rand - exact)
plt.title("error griddata")

plt.show()

【参考方案1】：

好问题！ （还有漂亮的情节！）

对于非结构化数据，您需要切换回用于非结构化数据的函数。 griddata 是一种选择，但使用三角测量和线性插值。这会导致三角形边界处的“硬”边缘。

样条是径向基函数。用 scipy 术语来说，你想要scipy.interpolate.Rbf。我建议在三次样条上使用function="linear" 或function="thin_plate"，但也可以使用三次。 （与线性或薄板样条相比，三次样条会加剧“过冲”问题。）

需要注意的是，径向基函数的这种特殊实现将始终使用数据集中的所有点。这是最准确和最平滑的方法，但随着输入观察点数量的增加，它的扩展性很差。有几种方法可以解决这个问题，但事情会变得更加复杂。我将把它留给另一个问题。

无论如何，这里有一个简化的例子。我们将生成随机数据，然后在规则网格上的点处对其进行插值。 （注意输入不在规则网格上，插值点也不需要。）

import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt
np.random.seed(1977)

x, y, z = np.random.random((3, 10))

interp = scipy.interpolate.Rbf(x, y, z, function='thin_plate')

yi, xi = np.mgrid[0:1:100j, 0:1:100j]

zi = interp(xi, yi)

plt.plot(x, y, 'ko')
plt.imshow(zi, extent=[0, 1, 1, 0], cmap='gist_earth')
plt.colorbar()

plt.show()

---

样条类型的选择

我选择"thin_plate" 作为样条线的类型。我们的输入观察点范围从 0 到 1（它们由 np.random.random 创建）。请注意，我们的内插值略高于 1 且远低于零。这就是“过冲”。

---

线性样条曲线将完全避免过冲，但您最终会得到“靶心”模式（尽管远没有 IDW 方法那么严重）。例如，这是用线性径向基函数插值的完全相同的数据。请注意，我们的插值永远不会高于 1 或低于 0：

---

高阶样条曲线会使数据中的趋势更加连续，但会出现更多过冲。默认的"multiquadric" 与薄板样条曲线非常相似，但会使事情变得更连续，并且过冲更糟：

---

但是，当您使用更高阶的样条曲线时，例如 "cubic"（三阶）：

和"quintic"（五阶）

只要稍微超出输入数据，您就真的会得到不合理的结果。

---

无论如何，这里有一个简单的例子来比较随机数据上的不同径向基函数：

import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt
np.random.seed(1977)

x, y, z = np.random.random((3, 10))
yi, xi = np.mgrid[0:1:100j, 0:1:100j]

interp_types = ['multiquadric', 'inverse', 'gaussian', 'linear', 'cubic',
                'quintic', 'thin_plate']

for kind in interp_types:
    interp = scipy.interpolate.Rbf(x, y, z, function=kind)
    zi = interp(xi, yi)

    fig, ax = plt.subplots()
    ax.plot(x, y, 'ko')
    im = ax.imshow(zi, extent=[0, 1, 1, 0], cmap='gist_earth')
    fig.colorbar(im)
    ax.set(title=kind)
    fig.savefig(kind + '.png', dpi=80)

plt.show()

【讨论】：

我已经用过RBF，效果非常好，建议至少尝试一下！

感谢您抽出宝贵时间发布详尽的答案！这正是我所需要的。好的，我要试试 Rbf。我很欣赏有关注意事项的信息，这将节省大量时间，而不必找出困难的方法。这将会非常好玩！对问题的“次要”部分有什么想法吗？

感谢@heltonbiker 的鼓励！

...我也喜欢你的脚本短小精悍。

当数据集中有很多点时，您能否详细说明如何使 RBF 更快？