使用Scipy/Numpy进行云点的二维插值，仅获取“有效”点

Question

使用Scipy/Numpy进行云点的二维插值，仅获取“有效”点

7

我有一个由人背部的摄影测量得到的云点。我正在尝试对其进行插值以获得常规网格，为此我使用了 scipy.interpolate ，目前效果良好。问题是：我正在使用的函数（scipy.interpolate.griddata）在平面x，y中使用云点的凸包，因此会产生一些不存在于原始曲面中的值，其周长具有凹面。

下图显示了原始云点（左侧）的示例（水平线所显示的内容实际上是点云密集的线状云），griddata 给我的结果在中间，而我想要获取的结果在右侧--即云点在x，y平面上的“阴影”，其中原始曲面中不存在的点将为零或Nan。

我知道可以从云点中移除Z坐标并检查每个网格位置的接近性，但这太过暴力，我相信这应该是点云应用程序上的常见问题。另一种可能性是在点云上执行一些numpy操作，找到一个numpy掩码或布尔2D数组以“应用”于来自griddata的结果，但我没有发现任何（这些操作超出了我的Numpy / Scipy知识范围）。

有什么建议吗？

谢谢您的阅读！

- heltonbiker

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- pv. · Accepted Answer

适当的口罩可以使用 KDTree 快速制造。griddata 使用的插值算法没有“有效”点的概念，因此您需要在插值之前或之后调整数据。

import numpy as np
from scipy.spatial import cKDTree as KDTree
from scipy.interpolate import griddata
import matplotlib.pyplot as plt

# Some input data
t = 1.2*np.pi*np.random.rand(3000)
r = 1 + np.random.rand(t.size)
x = r*np.cos(t)
y = r*np.sin(t)
z = x**2 - y**2

# -- Way 1: seed input with nan

def excluding_mesh(x, y, nx=30, ny=30):
    """
    Construct a grid of points, that are some distance away from points (x, 
    """

    dx = x.ptp() / nx
    dy = y.ptp() / ny

    xp, yp = np.mgrid[x.min()-2*dx:x.max()+2*dx:(nx+2)*1j,
                      y.min()-2*dy:y.max()+2*dy:(ny+2)*1j]
    xp = xp.ravel()
    yp = yp.ravel()

    # Use KDTree to answer the question: "which point of set (x,y) is the
    # nearest neighbors of those in (xp, yp)"
    tree = KDTree(np.c_[x, y])
    dist, j = tree.query(np.c_[xp, yp], k=1)

    # Select points sufficiently far away
    m = (dist > np.hypot(dx, dy))
    return xp[m], yp[m]

# Prepare fake data points
xp, yp = excluding_mesh(x, y, nx=35, ny=35)
zp = np.nan + np.zeros_like(xp)

# Grid the data plus fake data points
xi, yi = np.ogrid[-3:3:350j, -3:3:350j]
zi = griddata((np.r_[x,xp], np.r_[y,yp]), np.r_[z, zp], (xi, yi),
              method='linear')
plt.imshow(zi)
plt.show()

这个想法是在输入数据中添加包含nan值的虚假数据点。使用线性插值时，这些点将覆盖没有附近实际数据点的图像区域。

您也可以在插值后消除无效数据：

# -- Way 2: blot out afterward

xi, yi = np.mgrid[-3:3:350j, -3:3:350j]
zi = griddata((x, y), z, (xi, yi))

tree = KDTree(np.c_[x, y])
dist, _ = tree.query(np.c_[xi.ravel(), yi.ravel()], k=1)
dist = dist.reshape(xi.shape)
zi[dist > 0.1] = np.nan

plt.imshow(zi)
plt.show()