3

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

我有一个带有浮点值的 3D 直线网格（因此 3 个坐标轴 -1D numpy 数组 - 用于网格单元的中心和一个具有相应形状的 3D numpy 数组，每个单元中心都有一个值），我想插值（或您可以将其称为子采样）将整个数组转换为具有线性插值的子采样数组（例如，大小因子为 5）。到目前为止，我看到的所有方法都涉及 2D 和 1D 插值或 Id 宁愿不使用的 VTK 技巧（可移植性）。

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

提前感谢您的任何建议

编辑： 这是数据的样子，“d”是一个 3D 数组，表示 3D 单元格网格。每个单元格都有一个标量浮点值（在我的情况下为压力），“x”、“y”和“z”是三个 1D 数组，其中包含每个单元格的单元格的空间坐标（请参阅形状以及“x”数组的方式好像）

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 数据“插值”到那个大批。我不确定“插值”是否是正确的术语。下采样？平均？这是数据的二维切片：密度图

4

2 回答 2

6

这是使用scipy.interpolate.griddata在不规则网格上进行 3D 插值的示例。

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

在此处输入图像描述

于 2013-04-01T22:26:21.150 回答
1

简而言之：分别在每个维度上进行插值是正确的方法。

您可以简单地平均每个 5x5x5 立方体并返回结果。但是，如果您的数据应该是连续的，您应该明白这不是好的二次采样做法，因为它可能会导致混叠。（此外，您不能合理地将其称为“插值”！）

良好的重采样滤波器需要比重采样因子更宽，以避免混叠。由于您正在下采样，您还应该意识到您的重采样滤波器需要根据目标分辨率而不是原始分辨率进行缩放——为了正确插值，它可能需要是 5x5x5 的 4 或 5 倍立方体。这是很多样本 -20*20*20远远超过5*5*5......

因此，重采样的实际实现通常单独过滤每个维度的原因是它更有效。通过 3 次通过，您可以对每个输出样本使用更少的乘法/累加操作来评估您的滤波器。

于 2013-04-01T18:19:59.773 回答