给定一组在(X,Y,Z)坐标中定义的点,在任意(X,Y)处插值Z值

Question

给定(X,Y,Z)坐标中的一组点是表面上的点,我希望能够在任意(X,Y)坐标处插值Z值.我发现使用mlab.griddata在网格上插值可以获得一些成功,但我希望能够为任何(X,Y)坐标调用通用函数.

这组点形成一个大致半球形的表面.为了简化问题,我试图编写一种方法,在下面的x,y和z坐标定义的半球的已知点之间插值.虽然有一个分析解决方案可以找到完美球体的z = f(x,y),这样你就不需要插值,实际的点集将不是一个完美的球体,所以我们应该假设我们需要在未知(X,Y)坐标处插值.使用点数据链接到IPython笔记本

resolution = 10
u = np.linspace(-np.pi / 2, np.pi / 2, resolution)
v = np.linspace(0, np.pi, resolution)

U, V = np.meshgrid(u, v)

xs = np.sin(U) * np.cos(V) 
ys = np.sin(U) * np.sin(V)
zs = np.cos(U)

我一直在使用scipy.interpolate.interp2d,它"返回一个函数,其调用方法使用样条插值来查找新点的值."

def polar(xs, ys, zs, resolution=10):
    rs = np.sqrt(np.multiply(xs, xs) + np.multiply(ys, ys))
    ts = np.arctan2(ys, xs)
    func = interp2d(rs, ts, zs, kind='cubic')
    vectorized = np.vectorize(func)

    # Guesses
    ri = np.linspace(0, rs.max(), resolution)
    ti = np.linspace(0, np.pi * 2, resolution)

    R, T = np.meshgrid(ri, ti)
    Z = vectorized(R, T)
    return R * np.cos(T), R * np.sin(T), Z

不幸的是,我得到了非常奇怪的结果,类似于另一个尝试使用interp2d的 StackOverflow 用户.

到目前为止我发现的最大成功是使用反平方来估计(X,Y)处的Z值.但该函数对于估计Z = 0附近的Z值并不完美.

如果z = f(x, y)在(x,y,z)中给出一组点,我该怎么做？我在这里遗漏了什么......我是否需要更多的点云来可靠地估算表面上的值？

编辑:

这是我写完的功能.该函数采用输入数组xs, ys, zs和插值x, y使用scipy.interpolate.griddata,不需要常规网格.我确信有一种更聪明的方法可以做到这一点,并希望任何更新,但它的工作原理,我不关心性能.包括一个片段,以防将来帮助任何人.

def interpolate(x, y, xs, ys, zs):
    r = np.sqrt(x*x + y*y)
    t = np.arctan2(y, x)

    rs = np.sqrt(np.multiply(xs, xs) + np.multiply(ys, ys))
    ts = np.arctan2(ys, xs)

    rs = rs.ravel()
    ts = ts.ravel()
    zs = zs.ravel()

    ts = np.concatenate((ts - np.pi * 2, ts, ts + np.pi * 2))
    rs = np.concatenate((rs, rs, rs))
    zs = np.concatenate((zs, zs, zs))


    Z = scipy.interpolate.griddata((rs, ts), zs, (r, t))
    Z = Z.ravel()
    R, T = np.meshgrid(r, t)
    return Z

Answer 1

你是说你已经尝试过使用griddata. 那么为什么那不起作用呢？griddata如果新点的间隔不规则，也适用。例如，

# Definitions of xs, ys and zs
nx, ny = 20, 30
x = np.linspace(0, np.pi, nx)
y = np.linspace(0, 2*np.pi, ny)

X,Y = np.meshgrid(x, y)

xs = X.reshape((nx*ny, 1))
ys = Y.reshape((nx*ny, 1))

## Arbitrary definition of zs
zs = np.cos(3*xs/2.)+np.sin(5*ys/3.)**2

## new points where I want the interpolations
points = np.random.rand(1000, 2)

import scipy.interpolate
zs2 = scipy.interpolate.griddata(np.hstack((xs, ys)), zs, points)

这不是你所追求的吗？