"""
3次元周期データのFFTを用いた補間処理を提供します。
このモジュールは、フーリエ変換とゼロパディングを利用して、3次元空間に周期的に配置された離散データを高解像度化するための関数を含みます。
主に科学技術計算や画像処理の分野で、データの補間やアップサンプリングに利用されます。
:doc:`interpolate3d_fft_usage`
"""
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.fft import fftn, ifftn, fftshift, ifftshift
[ドキュメント]
def interpolate_3d_periodic_data_fft(data, interp_factor=(2, 2, 2)):
"""
3次元の等間隔周期データをFFT（高速フーリエ変換）で補間します。
この関数は、与えられた3次元データをフーリエ変換（FFT）により周波数領域に変換します。
周波数領域において、元のデータが占める領域以外の部分をゼロで埋める（ゼロパディング）ことで、
より高周波数成分がゼロであると仮定し、実空間での高解像度補間を実現します。
ゼロパディングされた周波数領域データに逆フーリエ変換（IFFT）を適用することで、
指定された補間倍率でサンプリングされた高解像度のデータが得られます。
補間は周期境界条件を仮定して行われます。
`fftshift` と `ifftshift` を使用して、FFTのDC成分が配列の中心にくるように処理することで、
ゼロパディングが適切に行われます。
:param data: np.ndarray
補間したい3次元の周期データ。形状は `(Nx, Ny, Nz)` である必要があります。
要素のデータ型は実数または複素数に対応します。
:param interp_factor: tuple
各次元 `(x, y, z)` の補間倍率を示す3つの整数 `(interp_x, interp_y, interp_z)` のタプル。
各要素は1以上の整数である必要があります。`1` を指定した場合、その次元は補間されません。
:returns: np.ndarray
補間された3次元データ。元のデータ型が実数であれば実数データが返されます。
この実装では最終的に `np.real` を適用しているため、実数部のみが返されます。
形状は `(Nx * interp_x, Ny * interp_y, Nz * interp_z)` となります。
:raises ValueError:
入力データが3次元のNumPy配列でない場合、または `interp_factor` が適切な形式でない場合に発生します。
"""
if not (isinstance(data, np.ndarray) and data.ndim == 3):
raise ValueError("入力データは3次元のNumPy配列である必要があります。")
if not (isinstance(interp_factor, tuple) and len(interp_factor) == 3 and all(isinstance(f, int) and f >= 1 for f in interp_factor)):
raise ValueError("interp_factorは各次元の整数補間倍率を示す3つの要素を持つタプルである必要があります。")
Nx, Ny, Nz = data.shape
interp_Nx, interp_Ny, interp_Nz = interp_factor
# 1. FFTを適用して周波数領域のデータを得る
F_data = fftshift(fftn(data))
# 2. 周波数領域でのゼロパディング
# 新しいサイズを計算
new_Nx = Nx * interp_Nx
new_Ny = Ny * interp_Ny
new_Nz = Nz * interp_Nz
# 新しいゼロパディングされた配列を初期化
F_data_padded = np.zeros((new_Nx, new_Ny, new_Nz), dtype=F_data.dtype)
# 元の周波数成分を新しい配列の中央にコピー
start_x = (new_Nx - Nx) // 2
end_x = start_x + Nx
start_y = (new_Ny - Ny) // 2
end_y = start_y + Ny
start_z = (new_Nz - Nz) // 2
end_z = start_z + Nz
F_data_padded[start_x:end_x, start_y:end_y, start_z:end_z] = F_data
# 3. 逆FFTを適用して補間されたデータを得る
interpolated_data = np.real(ifftn(ifftshift(F_data_padded))) * (new_Nx * new_Ny * new_Nz) / (Nx * Ny * Nz)
return interpolated_data
# --- 使用例 ---
if __name__ == "__main__":
# 1. サンプル3次元データの作成
Nx_orig, Ny_orig, Nz_orig = 10, 10, 10
x_orig = np.linspace(0, 2 * np.pi, Nx_orig, endpoint=False)
y_orig = np.linspace(0, 2 * np.pi, Ny_orig, endpoint=False)
z_orig = np.linspace(0, 2 * np.pi, Nz_orig, endpoint=False)
X_orig, Y_orig, Z_orig = np.meshgrid(x_orig, y_orig, z_orig, indexing='ij')
original_data = np.sin(X_orig * 2) + np.cos(Y_orig * 3) + np.sin(Z_orig * 1.5)
print(f"元のデータの形状: {original_data.shape}")
# 2. FFTによる補間
interp_factor = (4, 4, 4) # 各次元を4倍に補間
interpolated_data = interpolate_3d_periodic_data_fft(original_data, interp_factor)
print(f"補間後のデータの形状: {interpolated_data.shape}")
# **ここから修正点**
# interpolated_data の形状から新しい次元数を取得する
new_Nx, new_Ny, new_Nz = interpolated_data.shape
# **修正点ここまで**
# 3. 結果の可視化 (スライス表示)
slice_idx_orig_x = Nx_orig // 2
slice_idx_interp_x = new_Nx // 2 # 補間後のデータに対応するインデックス
fig = plt.figure(figsize=(14, 7))
ax1 = fig.add_subplot(121)
c1 = ax1.imshow(original_data[slice_idx_orig_x, :, :], origin='lower', cmap='viridis',
extent=[0, 2 * np.pi, 0, 2 * np.pi])
fig.colorbar(c1, ax=ax1, fraction=0.046, pad=0.04)
ax1.set_title(f'Original Data (X-slice at {slice_idx_orig_x})')
ax1.set_xlabel('Y-axis')
ax1.set_ylabel('Z-axis')
ax2 = fig.add_subplot(122)
c2 = ax2.imshow(interpolated_data[slice_idx_interp_x, :, :], origin='lower', cmap='viridis',
extent=[0, 2 * np.pi, 0, 2 * np.pi])
fig.colorbar(c2, ax=ax2, fraction=0.046, pad=0.04)
ax2.set_title(f'Interpolated Data (X-slice at {slice_idx_interp_x})')
ax2.set_xlabel('Y-axis')
ax2.set_ylabel('Z-axis')
plt.tight_layout()
plt.show()
# オプション：3D可視化（データ点が多すぎると描画が重くなる可能性があります）
fig_3d = plt.figure(figsize=(12, 6))
ax_orig_3d = fig_3d.add_subplot(121, projection='3d')
sc_orig = ax_orig_3d.scatter(X_orig.flatten(), Y_orig.flatten(), Z_orig.flatten(),
c=original_data.flatten(), cmap='viridis', s=20)
fig_3d.colorbar(sc_orig, ax=ax_orig_3d, shrink=0.5, aspect=5)
ax_orig_3d.set_title('Original Data Points')
ax_orig_3d.set_xlabel('X')
ax_orig_3d.set_ylabel('Y')
ax_orig_3d.set_zlabel('Z')
ax_interp_3d = fig_3d.add_subplot(122, projection='3d')
# **ここから修正点**
# new_Nx, new_Ny, new_Nz がここで定義されているため、これらを使用できます
x_interp = np.linspace(0, 2 * np.pi, new_Nx, endpoint=False)
y_interp = np.linspace(0, 2 * np.pi, new_Ny, endpoint=False)
# **修正点ここまで**
X_interp_slice, Y_interp_slice = np.meshgrid(x_interp, y_interp, indexing='ij')
z_slice_idx = new_Nz // 2 # 補間後のデータに対応するインデックス
surf_interp = ax_interp_3d.plot_surface(X_interp_slice, Y_interp_slice,
interpolated_data[:, :, z_slice_idx],
cmap='viridis', edgecolor='none')
fig_3d.colorbar(surf_interp, ax=ax_interp_3d, shrink=0.5, aspect=5)
ax_interp_3d.set_title(f'Interpolated Data (Z-slice at {z_slice_idx})')
ax_interp_3d.set_xlabel('X')
ax_interp_3d.set_ylabel('Y')
ax_interp_3d.set_zlabel('Value')
plt.tight_layout()
plt.show()