openpiv は何をしますか？
概要
OpenPIV（Open Particle Image Velocimetry）は、PIV画像ペアから流体の流れを解析します。前処理、相互相関、ベクトル検証、外れ値置換、平滑化、物理単位へのスケーリングをカバーしています。
以下はすべて openpiv 0.25.4 で検証されています。リリースごとにAPIが変わるため、異なるバージョンでスニペットを信頼する前に inspect.signature() を確認してください。
使用するタイミング
このスキルは、実験的なPIVや流れ可視化画像ペアを扱う場合に使用します：2次元速度場の測定、検査窓パラメータの調整、ベクトルの検証、渦度・ひずみ速度・乱流統計量の導出。流れを測定するのではなくシミュレーションする場合は、代わりにCFDスキルを使用してください。
クイックスタート
OpenPIVをインストール：
uv pip install openpiv
# 解析を再現可能にする必要がある場合は固定 -- 以下のすべてのスニペットはこのバージョンで確認されています。
uv pip install "openpiv==0.25.4"
画像ペアでPIV解析を実行：
import numpy as np
from openpiv import tools, pyprocess, validation, filters, scaling
frame_a = tools.imread("image_a.bmp")
frame_b = tools.imread("image_b.bmp")
# 相互相関。sig2noise_methodがNoneでない場合は(u, v, s2n)を返す。
u, v, s2n = pyprocess.extended_search_area_piv(
frame_a.astype(np.int32),
frame_b.astype(np.int32),
window_size=32,
overlap=12,
dt=0.02,
search_area_size=38,
correlation_method="linear", # search_area_size > window_sizeの場合に必要
sig2noise_method="peak2peak",
)
x, y = pyprocess.get_coordinates(
image_size=frame_a.shape,
search_area_size=38,
overlap=12,
)
# flagsはブール配列：Trueは異常ベクトルを示す。
flags = validation.sig2noise_val(s2n, threshold=1.05)
u, v = filters.replace_outliers(u, v, flags, method="localmean", max_iter=3, kernel_size=2)
# 物理単位にスケーリングし、プロットのために画像座標に反転。
x, y, u, v = scaling.uniform(x, y, u, v, scaling_factor=96.52)
x, y, u, v = tools.transform_coordinates(x, y, u, v)
tools.save("vectors.txt", x, y, u, v, flags)
または、バンドルされているCLIを使用します。これは上記のパイプラインをラップしています：
python skills/openpiv/scripts/runner.py \
--image frame_a.bmp --image frame_b.bmp --output_dir results --verbose
コアコンセプト
PIVの基礎
粒子画像流速測定法は、2枚の画像間で照明されたトレーサー粒子を追跡することにより流体速度を測定する光学的手法です。
処理フロー：
- 既知の時間間隔
dtで分離された画像ペア（frame_a、frame_b）を取得。 - 画像を検査窓に分割。
- 対応する窓を相互相関させ、ピーク変位を見つける。
- ベクトルを検証（信号対雑音比、グローバル範囲、局所メジアン）。
- 異常ベクトルを補間値で置換。
- ピクセル変位を物理単位にスケーリング。
検査窓パラメータ
window_size — 相関窓のサイズ（ピクセル、通常16～128）。窓が大きいほど相関は向上しますが、空間分解能は粗くなります。
overlap — 隣接する窓間で共有されるピクセル数（通常 window_size の50～75%）。オーバーラップが大きいとベクトル密度とコストが増加しますが、隣接ベクトルは独立ではなく相関を持つようになります。
search_area_size — 2枚目のフレームで探索される窓のサイズ。window_size 以上である必要があり、数ピクセル大きくすると大きな変位に対応できます。拡張探索領域は correlation_method="linear" と組み合わせて使用します — デフォルトの "circular" はFFTの折り返しに依存し、大きな変位を小さな変位にエイリアスします。references/advanced_algorithms.md を参照。
経験則：最大変位を window_size の約4分の1未満に保ち、窓あたり5～10個の粒子を目指します。
信号対雑音比
s2n は相関ピークがどの程度明確かを測定します。sig2noise_method はその計算方法を制御します — "peak2mean"（関数のデフォルト）または "peak2peak"。両者は異なるスケールであるため、一方に調整されたしきい値は他方では意味がありません。典型的な peak2peak のしきい値は1.05～1.3です。
flags = validation.sig2noise_val(s2n, threshold=1.05)
# flagsはbool：True == 異常。`~flags` で良ベクトルを選択。
一般的な操作
動的マスキング
マスキングは openpiv.preprocess にあり、openpiv.masking モジュールではありません。(image, mask) タプルを返し、float画像を期待します。
from openpiv import preprocess
# method="edges"：暗くて鋭いエッジのオブジェクト用；"intensity"：高コントラストのオブジェクト用。
frame_a_masked, mask_a = preprocess.dynamic_masking(
frame_a.astype(np.float64), method="intensity", filter_size=7, threshold=0.005
)
frame_b_masked, mask_b = preprocess.dynamic_masking(
frame_b.astype(np.float64), method="intensity", filter_size=7, threshold=0.005
)
返された画像を相関ステップに渡します — マスク領域は既にゼロになっています。元のフレームに mask を乗算しないでください：マスキングは既に適用されており、method="edges" の場合、マスクはブール値ではなく uint8 の0/255で返されるため、乗算すると画像が255倍にリスケールされます。
マルチパス処理
マルチパス（窓変形）は openpiv.windef にあり、PIVSettings データクラスによって駆動されます。pyprocess にはマルチパスのエントリポイントはありません。
import numpy as np
from openpiv import scaling, windef
settings = windef.PIVSettings()
settings.windowsizes = (64, 32, 16) # パスごとに1エントリ、減少（これがデフォルトでもある）
settings.overlap = (32, 16, 8) # windowsizesと同じ長さ
settings.num_iterations = 3 # 実際に実行するパス数
settings.sig2noise_threshold = 1.05
x, y, u, v, flags = windef.simple_multipass(
frame_a.astype(np.int32), frame_b.astype(np.int32), settings
)
# 出力はフレームあたりのピクセル数 — 自分で変換する。scaling.uniformは単にscaling_factorで割るだけなので、dtは別途適用する。
dt = 0.02
x, y, u, v = scaling.uniform(x, y, u, v, scaling_factor=96.52)
u, v = u / dt, v / dt
simple_multipass は既に検証、外れ値置換、残ったNaNをゼロで埋め、transform_coordinates を呼び出しているため、これらの手順を繰り返さないでください。
単位の落とし穴： PIVSettings には dt と scaling_factor フィールドがありますが、windef はどちらも使用しません — first_pass は dt なしで extended_search_area_piv を呼び出すため、マルチパスチェーン全体はフレームあたりのピクセル数で動作します。settings.dt = 0.02 と設定しても、返される値は変わりません。上記のように事後に変換します。
個々のパスを制御するには、windef.first_pass と windef.multipass_img_deform がより低レベルの構成要素です。
検証と後処理
検証手法
すべてのバリデータはブール配列を返し、Trueは異常ベクトルを示します。
# 信号対雑音
flags = validation.sig2noise_val(s2n, threshold=1.05)
# グローバル範囲 — (最小、最大) のタプルを位置引数または u_thresholds/v_thresholds として取る。
flags = validation.global_val(u, v, (-300, 300), (-300, 300))
# 局所メジアン — u_threshold と v_threshold は必須；size は近傍の半幅。
flags = validation.local_median_val(u, v, u_threshold=30.0, v_threshold=30.0, size=1)
# ブールORで結合（np.maximumではない — これらはブール配列）。
flags = (
validation.sig2noise_val(s2n, threshold=1.05)
| validation.global_val(u, v, (-300, 300), (-300, 300))
| validation.local_median_val(u, v, u_threshold=30.0, v_threshold=30.0)
)
これらのしきい値は u と v の単位で設定し、フレームあたりのピクセル数では設定しないでください。
extended_search_area_piv は dt で除算するため、dt=0.02 の場合、3 px/フレームの変位は150 px/秒になります。上記のしきい値はそのケースに適しています；PIV文献や PIVSettings.min_max_u_disp で使用される (-30, 30) の数値はpx/フレームの制限であり、これをpx/秒の出力に適用するとフィールド全体が拒否されます。スケーリング前に検証するか、しきい値も 1/dt でスケーリングします。
外れ値置換
u, v = filters.replace_outliers(
u, v, flags, method="localmean", max_iter=3, tol=1e-3, kernel_size=2
)
method は "localmean"、"disk"、"distance" の3つだけを受け入れます。認識されない名前は拒否されず、すべてゼロのカーネルにフォールスルーし、役に立たないフィールドを静かに返します。置換はフラグが付けられた位置を補間値で埋めることに注意してください — その後NaNで上書きすると、置換は無駄になります。どちらか一方を選択します：
# 補間する代わりに、フラグ付きベクトルを解析から完全に除外する。
u = np.where(flags, np.nan, u)
v = np.where(flags, np.nan, v)
平滑化
平滑化は openpiv.smoothn.smoothn です。openpiv.smooth モジュールはありません。最初の要素が平滑化されたフィールドであるタプルを返し、NaN入力は受け入れません。
from openpiv.smoothn import smoothn
u_smooth, *_ = smoothn(np.nan_to_num(u), s=0.5) # s：大きいほど滑らか
v_smooth, *_ = smoothn(np.nan_to_num(v), s=0.5)
u_smooth = np.asarray(u_smooth)
可視化
ベクトル場プロット
display_vector_field は保存されたベクトルファイルを読み取り、内部で plt.show() を呼び出すため、バッチ実行には非対話型バックエンドを選択します。
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from openpiv import tools
fig, ax = plt.subplots(figsize=(8, 8))
tools.display_vector_field(
"vectors.txt",
ax=ax,
scaling_factor=96.52, # scaling.uniformで使用したものと同じ係数、画像にマッピングし直すため
scale=50,
width=0.0035,
on_img=True,
image_name="frame_a.bmp",
)
fig.savefig("vector_field.png", dpi=150, bbox_inches="tight")
plt.close(fig)
カスタム可視化
import numpy as np
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
mag = np.sqrt(u**2 + v**2)
for ax, field, title, cmap in [
(axes[0], mag, "Velocity Magnitude", "viridis"),
(axes[1], u, "U Velocity", "RdBu_r"),
(axes[2], v, "V Velocity", "RdBu_r"),
]:
im = ax.imshow(field, cmap=cmap)
ax.set_title(title)
plt.colorbar(im, ax=ax)
fig.tight_layout()
fig.savefig("velocity_components.png")
plt.close(fig)
解析関数
scripts/analyze.py はこれらを runner.py によって書き込まれた params.npz に対してバンドルしています。保存された座標から物理的なグリッド間隔を推測するため、微分は単位長さあたりで出力されます：
import sys
sys.path.insert(0, "skills/openpiv/scripts")
from analyze import PIVAnalyzer
piv = PIVAnalyzer("results/params.npz")
vorticity = piv.compute_vorticity() # dv/dx - du/dy
exx, eyy, exy = piv.compute_strain()
stats = piv.compute_statistics() # u_mean, v_mean, rms_u, rms_v, tke
piv.plot_vector_field(save_path="quiver.png")
インラインで計算したい場合のスタンドアロン形式：
渦度
def compute_vorticity(u, v, dx=1.0, dy=None):
"""面外渦度 dv/dx - du/dy。物理的なグリッド間隔を渡し、1.0ではない。"""
dy = dx if dy is None else dy
return np.gradient(v, dx, axis=1) - np.gradient(u, dy, axis=0)
グリッド間隔は物理単位で (window_size - overlap) / scaling_factor であるため、dx=1.0 のままにすると、単位長さあたりではなくグリッドセルあたりの渦度が得られます。
符号の慣例： runner.py は transform_coordinates で終了し、グリッドを右手系のy上向きフレームにラベル付けし直しますが、行は画像の順序のままであるため、保存された y は行インデックスが大きくなるにつれて減少します。上記のスタンドアロン形式はその逆を仮定しているため、params.npz フィールドでは -du/dy を返し、渦度とせん断ひずみの符号が反転します — axis=0 の微分を否定するか、保存された座標から方向を読み取る PIVAnalyzer を使用します。
ひずみ速度
def compute_strain(u, v, dx=1.0, dy=None):
"""2Dひずみ速度テンソルの (exx, eyy, exy) を返す。"""
dy = dx if dy is None else dy
du_dx = np.gradient(u, dx, axis=1)
du_dy = np.gradient(u, dy, axis=0)
dv_dx = np.gradient(v, dx, axis=1)
dv_dy = np.gradient(v, dy, axis=0)
return du_dx, dv_dy, 0.5 * (du_dy + dv_dx)
乱流統計量
def compute_statistics(u, v):
"""単一フレームの空間統計量。レイノルズ分解ではない。"""
u_prime = u - np.nanmean(u)
v_prime = v - np.nanmean(v)
rms_u, rms_v = np.nanstd(u_prime), np.nanstd(v_prime)
return {
"u_mean": np.nanmean(u),
"v_mean": np.nanmean(v),
"rms_u": rms_u,
"rms_v": rms_v,
"tke": 0.5 * (rms_u**2 + rms_v**2),
}
注意： 1フレームの空間平均を減算すると空間分散が測定され、これは均質なフィールドに対してのみ乱流強度と等しくなります。真のレイノルズ分解には画像ペアのアンサンブルが必要です：時間軸で平均し、その平均場を各実現から差し引きます。
CLIの使い方
# 基本実行
python skills/openpiv/scripts/runner.py \
--image img1.bmp --image img2.bmp --output_dir results --verbose
# 動的マスキングによるチューニングパラメータ
python skills/openpiv/scripts/runner.py \
--image frame_a.bmp \
--image frame_b.bmp \
--output_dir results \
--window_size 32 \
--overlap 12 \
--search_area 38 \
--dt 0.02 \
--scaling 96.52 \
--threshold 1.05 \
--mask dynamic \
--mask_method intensity \
--verbose
CLIオプション
| オプション | デフォルト | 説明 |
|---|---|---|
--image | 必須 | 画像ファイル；ペアのため正確に2回指定 |
--output_dir | results | 出力ディレクトリ（存在しない場合は作成） |
--window_size | 32 | 検査窓サイズ（px） |
--overlap | 12 | 窓のオーバーラップ（px） |
--search_area | 38 | 探索領域サイズ（px）。--window_size 以上である必要がある |
--dt | 0.02 | フレーム間の時間（秒） |
--scaling | 96.52 | スケーリング係数、物理単位あたりのピクセル数（例：px/mm） |
--threshold | 1.05 | peak2peak 信号対雑音しきい値 |
--mask | none | none または dynamic（openpiv.preprocess.dynamic_masking） |
--mask_method | intensity | edges または intensity。--mask dynamic の場合に使用 |