3

私がやりたいことはかなり単純ですが、これまでのところ簡単なアプローチは見つかりませんでした:

float 値を持つ 3D 直線グリッド (したがって、グリッド セルの中心には 3 つの座標軸 - 1D numpy 配列 - と、各セルの中心に値を持つ対応する形状を持つ 3D numpy 配列) があり、補間したい (またはサブサンプルと呼ぶことができます) この配列全体を、線形補間を使用してサブサンプリングされた配列 (たとえば、サイズ係数 5) に変換します。私がこれまで見てきたすべてのアプローチには、2D と 1D 補間または使用したくない VTK トリック (移植性) が含まれます。

3D配列で同時に5x5x5セルを取得し、平均して各方向に5倍小さい配列を返すのと同等のアプローチを誰かが提案できますか?

提案をお寄せいただきありがとうございます

編集: データは次のようになります。「d」は、セルの 3D グリッドを表す 3D 配列です。各セルにはスカラー浮動小数点値 (私の場合は圧力) があり、「x」、「y」、および「z」は、すべてのセルのセルの空間座標を含む 3 つの 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 分の 1 に縮小) で 3D 配列を作成することです。配列。ただし、「補間」が正しい用語であるかどうかはわかりません。ダウンサンプリング？平均化？以下は、データの 2D スライスです。密度マップ

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*205*5*5

したがって、リサンプリングの実際の実装が通常、各次元を個別にフィルタリングする理由は、その方が効率的だからです。3 つのパスを取得することで、出力サンプルあたりの乗算/累算操作をはるかに少なくして、フィルターを評価できます。

于 2013-04-01T18:19:59.773 に答える