仅获得“有效"的使用Scipy/Numpy在浊点的2D插值中的点 [英] Get only "valid" points in 2D interpolation of cloud point using Scipy/Numpy

问题描述

我有一个通过摄影测量法从一个人的背上得到的浊点.我正在尝试对其进行插值以获得规则的网格，为此，到目前为止，我正在使用 scipy.interpolate 并获得良好的结果.问题是:我正在使用的函数( scipy.interpolate.griddata )使用平面x，y上的浊点的凸包，因此得出的一些值不存在于原始表面具有凹入的周长.

I have a cloud point obtained from photogrammetry from a person's back. I'm trying to interpolate it to get a regular grid, and for that I'm using scipy.interpolate with good results so far. The problem is: the function I'm using (scipy.interpolate.griddata) uses the convex hull of the cloudpoint in the plane x,y, thus giving as result some values that don't exist in the original surface, which has a concave perimeter.

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

The following illustration shows the original cloudpoint at the left (what is displayed as horizontal lines is actually a dense line-shaped cloud of points), the result that griddata gives me in the middle, and the result I would like to get at the right -- kind of the "shadow" of the cloudpoint on the x,y plane, where non-existing points in the original surface would be zeros or Nans.

我知道我可以删除浊点上的Z坐标并检查每个网格位置是否接近，但这是蛮力的，我相信这应该是点云应用程序中的常见问题.另一种可能是在点云上执行一些numpy操作，找到一个numpy掩码或布尔2D数组以应用"来自 griddata 的结果，但是我没有找到任何(这些操作有点超出了我的Numpy/Scipy知识).

I know I could remove the Z coordinate on the cloudpoint and check each grid position for proximity, but this is so brute-force, and I believe this should be a common problem on point-cloud applications. Another possibility could be some numpy operation to perform over the point-cloud, finding a numpy mask or boolean 2D-array to "apply" over the result from griddata, but I didn't find any (these operations are a bit beyond my Numpy/Scipy knowledge).

有什么建议吗?

感谢阅读！

推荐答案

可以使用 KDTree 快速构建合适的蒙版.griddata使用的插值算法没有有效"点的概念，因此您需要在插值之前或之后调整数据.

Suitable masks can be constructed fast using KDTree. The interpolation algorithm used by griddata does not have a notion of "valid" points, so you need to adjust your data either before or after the interpolation.

之前:

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 值的假数据点播种"输入数据.使用线性插值时，这些将遮挡图像附近没有实际数据点的区域.

The idea is to "seed" the input data with fake data points containing nan values. When using linear interpolation, these will blot out areas of the image that have no actual data points nearby.

您还可以在插值后清除无效数据:

You can also blot out invalid data after the interpolation:

# -- 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()

这篇关于仅获得“有效"的使用Scipy/Numpy在浊点的2D插值中的点的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！