Obtener solo puntos "válidos" en la interpolación 2D del punto de nube usando Scipy/Numpy

Question

Tengo un punto de nube obtenido de la fotogrametría desde la espalda de una persona. Estoy tratando de interpolarlo para obtener una grilla regular, y para eso estoy usando scipy.interpolate con buenos resultados hasta ahora. El problema es que la función que estoy usando (scipy.interpolate.griddata) utiliza el casco convexo del punto de nube en el plano x, y, dando como resultado algunos valores que no existen en la superficie original, que tiene un perímetro cóncavo.

La siguiente ilustración muestra el punto de niebla original a la izquierda (lo que se muestra como líneas horizontales es en realidad una nube de puntos en forma de línea densa), el resultado que griddata me da en el medio, y el resultado me gustaría llegar a la derecha - tipo de la "sombra" del punto de nube en el plano x, y, donde los puntos no existentes en la superficie original serían cero o Nans.

Sé que podría eliminar la coordenada Z del punto de turbidez y comprobar cada posición en la parrilla de proximidad, pero esto es por lo que la fuerza bruta, y creo que esto debería ser un problema común en las aplicaciones de nubes de puntos . Otra posibilidad podría ser alguna operación numpy para realizar sobre la nube de puntos, encontrar una máscara numpy o boolean 2D-array para "aplicar" sobre el resultado de griddata, pero no encontré ninguna (estas operaciones son un poco más allá de mi Numpy/Scipy conocimiento).

¿Alguna sugerencia?

Gracias por leer!

Answer 1

Las máscaras adecuadas se pueden construir rápidamente usando KDTree. El algoritmo de interpolación utilizado por griddata no tiene una noción de puntos "válidos", por lo que debe ajustar sus datos antes o después de la interpolación.

Antes:

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

la idea es "semilla" los datos de entrada con puntos de datos falsos que contienen nan valores. Al usar la interpolación lineal, estos borrarán las áreas de la imagen que no tienen puntos de datos reales cerca.

También puede borrar datos no válidos después de la interpolación:

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