给定在 (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) 坐标处插入值。 Link to IPython notebook with point data

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

不幸的是，我得到了非常奇怪的结果，类似于另一个 StackOverflow user who tried to use interp2d。

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

如果给定 (x, y, z) 中的一组点，我该怎么做才能获得函数 z = f(x, y)？我在这里遗漏了什么...我需要的不仅仅是点云来可靠地估计表面上的值吗？

编辑：

这是我最终编写的函数。该函数采用xs, ys, zs 的输入数组并使用scipy.interpolate.griddata 在x, y 处进行插值，这不需要常规网格。我确信有一种更聪明的方法可以做到这一点，并且会感谢任何更新，但它有效，我不关心性能。包括一个 sn-p，以防将来对任何人有所帮助。

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)

这不是你想要的吗？

Answer 2

如果我理解您的问题，您的点 xs、ys、zs 定义为

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

您想要的是能够对给定的 x 和 y 进行插值并找到 z 值？为什么需要插值？上面的方程代表一个球体，可以改写为xs*xs + ys*ys + zs*zs = 1，所以这个问题有一个简单的解析解：

def Z(X, Y):
    return np.sqrt(1-X**2-Y**2)
    ## or return -np.sqrt(1-X**2-Y**2) since this equation has two solutions

除非我误解了这个问题。