Interpolación / submuestreo de datos 3D en python sin VTK

Lo que quiero hacer es bastante simple, pero hasta ahora no he encontrado un enfoque directo:

Tengo una cuadrícula rectilínea 3D con valores flotantes (por lo tanto, 3 ejes de coordenadas -1D matrices numpy- para los centros de las celdas de la cuadrícula y una matriz numpy 3D con la forma correspondiente con un valor para cada centro de celda) y quiero interpolar (o puede llamarla submuestra) esta matriz completa a una matriz submuestreada (por ejemplo, factor de tamaño de 5) con interpolación lineal. Todos los enfoques que he visto hasta ahora implican interpolación en 2D y luego en 1D o trucos VTK que Id no prefiere usar (portabilidad).

¿Podría alguien sugerir un enfoque que sería equivalente a tomar 5x5x5 celdas al mismo tiempo en la matriz 3D, promediando y devolviendo una matriz 5 veces más pequeña en cada dirección?

Gracias de antemano por cualquier sugerencia

EDITAR: Así es como se ven los datos, ‘d’ es una matriz 3D que representa una cuadrícula 3D de celdas. Cada celda tiene un valor flotante escalar (presión en mi caso) y ‘x’, ‘y’ y ‘z’ son tres matrices 1D que contienen las coordenadas espaciales de las celdas de cada celda (vea las formas y cómo la matriz ‘x’ parece)

    In [42]: x.shape Out[42]: (181L,) In [43]: y.shape Out[43]: (181L,) In [44]: z.shape Out[44]: (421L,) In [45]: d.shape Out[45]: (181L, 181L, 421L) In [46]: x Out[46]: array([-0.410607 , -0.3927568 , -0.37780656, -0.36527296, -0.35475321, -0.34591168, -0.33846866, -0.33219107, -0.32688467, -0.3223876 , ... 0.34591168, 0.35475321, 0.36527296, 0.37780656, 0.3927568 , 0.410607 ]) 

    Lo que quiero hacer es crear una matriz 3D con, digamos, una forma de 90x90x210 (aproximadamente una reducción de tamaño por un factor de 2) al sub-muestrear primero las coordenadas de los ejes en las matrices con las dimensiones anteriores y luego “interpolar” los datos 3D para que formación. No estoy seguro si “interpolar” es el término correcto sin embargo. Downsampling? Averaging? Aquí hay una porción 2D de los datos: mapas de densidad

    Aquí hay un ejemplo de interpolación 3D en una cuadrícula irregular utilizando scipy.interpolate.griddata .

     import numpy as np import scipy.interpolate as interpolate import matplotlib.pyplot as plt def func(x, y, z): return x ** 2 + y ** 2 + z ** 2 # Nx, Ny, Nz = 181, 181, 421 Nx, Ny, Nz = 18, 18, 42 subsample = 2 Mx, My, Mz = Nx // subsample, Ny // subsample, Nz // subsample # Define irregularly spaced arrays x = np.random.random(Nx) y = np.random.random(Ny) z = np.random.random(Nz) # Compute the matrix D of shape (Nx, Ny, Nz). # D could be experimental data, but here I'll define it using func # D[i,j,k] is associated with location (x[i], y[j], z[k]) X_irregular, Y_irregular, Z_irregular = ( x[:, None, None], y[None, :, None], z[None, None, :]) D = func(X_irregular, Y_irregular, Z_irregular) # Create a uniformly spaced grid xi = np.linspace(x.min(), x.max(), Mx) yi = np.linspace(y.min(), y.max(), My) zi = np.linspace(y.min(), y.max(), Mz) X_uniform, Y_uniform, Z_uniform = ( xi[:, None, None], yi[None, :, None], zi[None, None, :]) # To use griddata, I need 1D-arrays for x, y, z of length # len(D.ravel()) = Nx*Ny*Nz. # To do this, I broadcast up my *_irregular arrays to each be # of shape (Nx, Ny, Nz) # and then use ravel() to make them 1D-arrays X_irregular, Y_irregular, Z_irregular = np.broadcast_arrays( X_irregular, Y_irregular, Z_irregular) D_interpolated = interpolate.griddata( (X_irregular.ravel(), Y_irregular.ravel(), Z_irregular.ravel()), D.ravel(), (X_uniform, Y_uniform, Z_uniform), method='linear') print(D_interpolated.shape) # (90, 90, 210) # Make plots fig, ax = plt.subplots(2) # Choose az value in the uniform z-grid # Let's take the middle value zindex = Mz // 2 z_crosssection = zi[zindex] # Plot a cross-section of the raw irregularly spaced data X_irr, Y_irr = np.meshgrid(sorted(x), sorted(y)) # find the value in the irregular z-grid closest to z_crosssection z_near_cross = z[(np.abs(z - z_crosssection)).argmin()] ax[0].contourf(X_irr, Y_irr, func(X_irr, Y_irr, z_near_cross)) ax[0].scatter(X_irr, Y_irr, c='white', s=20) ax[0].set_title('Cross-section of irregular data') ax[0].set_xlim(x.min(), x.max()) ax[0].set_ylim(y.min(), y.max()) # Plot a cross-section of the Interpolated uniformly spaced data X_unif, Y_unif = np.meshgrid(xi, yi) ax[1].contourf(X_unif, Y_unif, D_interpolated[:, :, zindex]) ax[1].scatter(X_unif, Y_unif, c='white', s=20) ax[1].set_title('Cross-section of downsampled and interpolated data') ax[1].set_xlim(x.min(), x.max()) ax[1].set_ylim(y.min(), y.max()) plt.show() 

    introduzca la descripción de la imagen aquí

    En resumen: hacer la interpolación en cada dimensión por separado es el camino correcto.


    Simplemente puede promediar cada cubo de 5x5x5 y devolver los resultados. Sin embargo, si se supone que sus datos son continuos, debe comprender que no es una buena práctica de submuestreo, ya que es probable que provoque aliasing. (Además, ¡no puedes llamarlo “interpolación”!)

    Los buenos filtros de remuestreo deben ser más anchos que el factor de remuestreo para evitar el alias. Como está reduciendo el muestreo, también debe darse cuenta de que su filtro de remuestreo debe escalarse de acuerdo con la resolución de destino , no con la resolución original. Para interpolar adecuadamente, es probable que tenga que ser 4 o 5 veces más ancho que su 5x5x5 cubo. Esta es una gran cantidad de muestras: 20*20*20 es mucho más que 5*5*5

    Entonces, la razón por la que las implementaciones prácticas de remuestreo típicamente filtran cada dimensión por separado es que es más eficiente. Al tomar 3 pases, puede evaluar su filtro utilizando muchas menos operaciones de multiplicación / acumulación por muestra de salida.