没有VTK的python中3D数据的插值/子采样 [英] Interpolation/subsampling of 3D data in python without VTK

问题描述

我想做的很简单，但到目前为止我还没有找到一个简单的方法:

What I want to do is rather simple but I havent found a straightforward approach thus far:

我有一个带有浮点值的 3D 直线网格(因此有 3 个坐标轴 -1D numpy 数组 - 用于网格单元的中心和一个具有相应形状的 3D numpy 数组，每个单元中心都有一个值)，我想使用线性插值将整个数组插入(或者您可以称之为子采样)到子采样数组(例如大小因子为 5).到目前为止，我所看到的所有方法都涉及 2D，然后是 1D 插值或 VTK 技巧，我不想使用这些技巧(可移植性).

I have a 3D rectilinear grid with float values (therefore 3 coordinate axes -1D numpy arrays- for the centers of the grid cells and a 3D numpy array with the corresponding shape with a value for each cell center) and I want to interpolate (or you may call it subsample) this entire array to a subsampled array (e.g. size factor of 5) with linear interpolation. All the approaches I've seen this far involve 2D and then 1D interpolation or VTK tricks which Id rather not use (portability).

有人可以提出一种方法，相当于在 3D 阵列中同时获取 5x5x5 个单元，平均并返回一个在每个方向上小 5 倍的阵列吗?

Could someone suggest an approach that would be the equivalent of taking 5x5x5 cells at the same time in the 3D array, averaging and returning an array 5times smaller in each direction?

提前感谢您的任何建议

下面是数据的样子，d"是一个 3D 数组，表示一个 3D 单元格网格.每个单元格都有一个标量浮点值(在我的例子中是压力)，'x'、'y' 和 'z' 是三个一维数组，包含每个单元格的单元格的空间坐标(参见形状和 'x' 数组如何看起来像)

Here's what the data looks like, 'd' is a 3D array representing a 3D grid of cells. Each cell has a scalar float value (pressure in my case) and 'x','y' and 'z' are three 1D arrays containing the spatial coordinates of the cells of every cell (see the shapes and how the 'x' array looks like)

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

我想要做的是创建一个 3D 数组，可以说形状为 90x90x210(大约缩小 2 倍)，方法是首先对具有上述维度的数组上的轴坐标进行二次采样，然后插值"3D数据到该数组.我不确定插值"是否是正确的术语.降采样?平均?这是数据的 2D 切片:

What I want to do is create a 3D array with lets say a shape of 90x90x210 (roughly downsize by a factor of 2) by first subsampling the coordinates from the axes on arrays with the above dimensions and then 'interpolating' the 3D data to that array. Im not sure whether 'interpolating' is the right term though. Downsampling? Averaging? Here's an 2D slice of the data:

推荐答案

以下是使用 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 a z 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()

这篇关于没有VTK的python中3D数据的插值/子采样的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！