只获取“有效"使用 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 掩码或布尔二维数组来应用"来自 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屋！