3d - 使用 Scipy/Numpy 在云点的 2D 插值中仅获取 "valid"个点

Question

我从一个人的背部通过摄影测量获得了一个浊点。我正在尝试对其进行插值以获得常规网格，为此我正在使用 scipy.interpolate到目前为止效果很好。问题是:我正在使用的函数( scipy.interpolate.griddata )使用平面 x,y 中云点的凸包，因此给出了一些原始表面中不存在的值，该表面具有凹周长.

下图显示了左边的原始浊点(显示为水平线的实际上是密集的线状点云)，结果griddata给我在中间，我想得到的结果是正确的 - 一种 x,y 平面上云点的“阴影”，其中原始表面中不存在的点将为零或 Nans。



我知道我可以删除云点上的 Z 坐标并检查每个网格位置的接近度，但这太暴力了，我相信这应该是点云应用程序的常见问题。
另一种可能性可能是在点云上执行一些 numpy 操作，找到一个 numpy 掩码或 bool 二维数组来“应用”来自 griddata 的结果。 ，但我没有找到任何(这些操作有点超出我的 Numpy/Scipy 知识)。

有什么建议吗？

谢谢阅读!

Answer 1

可以使用 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()