Scipy/Numpyを使用して雲点の2D補間で "有効な"点のみを取得します

Question

人の背中からの写真測量から得られた曇り点があります。私は通常のグリッドを取得するためにそれを補間しようとしています。そのために私はscipy.interpolateを使用しています。問題は、私が使用している関数（scipy.interpolate.griddata）は、x、y平面上のクラウドポイントの凸包を使用しているため、元のサーフェスには存在しないいくつかの値が得られます。

次の図は、水平線として表示されている元の曇り点（実際には濃い線状の雲の点です）を示しています。結果はgriddataです。オリジナルのサーフェスに存在しない点がゼロまたはナンズになるx、y平面上の曇り点の「影」の右の種類を取得します。

私はZは曇り点上の座標と近接のために各グリッドの位置を確認し削除することもできます知っているが、これはそうブルートフォースであり、私はこれがポイントクラウドアプリケーションに共通の問題であるべきと考えています。 もう1つの可能性は、ポイントクラウド上で実行するnumpy操作で、num33マスクまたはブール値の2次元配列を見つけてgriddataの結果に適用することができますが、何も見つかりませんでした。ナンシー/シフィー知識）。

提案がありますか？

読んでいただきありがとうございます！

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