Получите только «действительные» точки в 2D-интерполяции точки облачности с помощью Scipy/Numpy

Question

У меня есть точка облачности, полученная от фотограмметрии от спины человека. Я пытаюсь интерполировать его, чтобы получить регулярную сетку, и для этого я использую scipy.interpolate с хорошими результатами. Проблема в том, что функция, которую я использую (scipy.interpolate.griddata), использует выпуклую оболочку точки облачности в плоскости x, y, что дает в результате некоторые значения, которые не существуют на исходной поверхности, которая имеет вогнутый периметр.

На следующем рисунке показана исходная точка облачности слева (то, что отображается в виде горизонтальных линий, на самом деле представляет собой плотное линейное облако точек), результат griddata дает мне посередине, и результат, который я хотел бы получить справа - вид «тени» точки облачности на плоскости x, y, где несуществующие точки на исходной поверхности будут нулями или Nans.

Я знаю, что я мог бы удалить Z координату на cloudpoint и проверить каждую позицию сетки для близости, но это так перебор, и я считаю, что это должно быть общей проблема на точки облачных приложений , Другой возможностью может быть некоторая операция numpy для выполнения над облаком точек, нахождение маски numpy или логического 2D-массива для «применения» по результату от 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()