openpiv

OpenPIV

概要

OpenPIV（Open Particle Image Velocimetry）は、PIV 画像ペアから流体の流れを解析します。前処理、相互相関、ベクトルの検証、外れ値の置換、平滑化、物理単位へのスケーリングに対応しています。

以下の内容はすべて openpiv 0.25.4 に照らして検証済みです。API はリリース間で変更されるため、異なるバージョンでスニペットを使用する前に inspect.signature() を確認してください。

使用する場面

実験用 PIV または流れの可視化画像ペアを扱う場合に、このスキルを使用してください。2D 速度場の測定、検査窓パラメーターの調整、ベクトルの検証、渦度、ひずみ速度、乱流統計量の導出などに使えます。流れの測定ではなく シミュレーション を行う場合は、代わりに CFD スキルを使用してください。

クイックスタート

OpenPIV をインストールします。

uv pip install openpiv

# Pin it when the analysis needs to be reproducible -- this is the version every
# snippet below was checked against.
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")

# Cross-correlate. Returns (u, v, s2n) whenever sig2noise_method is not None.
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",   # required for 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 is a boolean array: True marks a spurious vector.
flags = validation.sig2noise_val(s2n, threshold=1.05)
u, v = filters.replace_outliers(u, v, flags, method="localmean", max_iter=3, kernel_size=2)

# Scale to physical units, then flip to image coordinates for plotting.
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 の基礎

粒子画像流速測定法（Particle Image Velocimetry）は、2 枚の画像間で照明されたトレーサー粒子を追跡し、流体速度を測定する光学的手法です。

処理の流れ：

1. 既知の時間間隔 dt を空けて、画像ペア（frame_a、frame_b）を撮影します。
2. 画像を検査窓に分割します。
3. 対応する窓を相互相関し、ピーク変位を求めます。
4. ベクトルを検証します（信号対雑音比、全体範囲、局所中央値）。
5. 異常なベクトルを補間値で置き換えます。
6. ピクセル変位を物理単位にスケーリングします。

検査窓のパラメーター

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 未満に抑え、1 窓あたり 5～10 個の粒子を目安にしてください。

信号対雑音比

s2n は、相関ピークの明瞭さを測定します。sig2noise_method はその計算方法を制御します。"peak2mean"（関数のデフォルト）または "peak2peak" を指定します。この 2 つは尺度が異なるため、一方に合わせて調整したしきい値は、もう一方では意味を成しません。peak2peak の一般的なしきい値は 1.05～1.3 です。

flags = validation.sig2noise_val(s2n, threshold=1.05)
# flags is bool: True == spurious. `~flags` selects the good vectors.

よくある操作

動的マスキング

マスキングは openpiv.masking モジュールではなく、openpiv.preprocess にあります。(image, mask) タプルを返し、浮動小数点画像を想定しています。

from openpiv import preprocess

# method="edges" for dark, sharp-edged objects; "intensity" for high-contrast objects.
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 倍にスケーリングされます。

マルチパス処理

マルチパス（窓変形）は PIVSettings データクラスを使って openpiv.windef に実装されています。pyprocess にマルチパス用のエントリーポイントはありません。

import numpy as np
from openpiv import scaling, windef

settings = windef.PIVSettings()
settings.windowsizes = (64, 32, 16)   # one entry per pass, decreasing (this is also the default)
settings.overlap = (32, 16, 8)        # same length as windowsizes
settings.num_iterations = 3           # number of passes to actually run
settings.sig2noise_threshold = 1.05

x, y, u, v, flags = windef.simple_multipass(
    frame_a.astype(np.int32), frame_b.astype(np.int32), settings
)

# Output is in PIXELS PER FRAME -- convert yourself. scaling.uniform only divides
# by scaling_factor, so apply dt separately.
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 は外れベクトルを示します。

# Signal-to-noise
flags = validation.sig2noise_val(s2n, threshold=1.05)

# Global range -- takes (min, max) TUPLES, positionally or as u_thresholds/v_thresholds.
flags = validation.global_val(u, v, (-300, 300), (-300, 300))

# Local median -- u_threshold and v_threshold are REQUIRED; size is the neighbourhood half-width.
flags = validation.local_median_val(u, v, u_threshold=30.0, v_threshold=30.0, size=1)

# Combine with boolean OR (not np.maximum -- these are bool arrays).
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/frame の変位は 150 px/s として渡されます。上記のしきい値はこのケースに適しています。PIV の文献や PIVSettings.min_max_u_disp で使われている (-30, 30) は px/frame の上限であり、 px/s の出力に適用すると、フィールド全体が除外されます。スケーリング前に検証するか、 しきい値も 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 で上書きすると、置換が 無駄になります。どちらか一方を選択してください。

# Keep flagged vectors out of the analysis entirely, instead of interpolating them.
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: larger == smoother
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,   # same factor used in scaling.uniform, to map back onto the image
    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):
    """Out-of-plane vorticity dv/dx - du/dy. Pass the physical grid spacing, not 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):
    """Return (exx, eyy, exy) of the 2D strain-rate tensor."""
    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):
    """Single-frame spatial statistics. NOT Reynolds decomposition."""
    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 の使用方法

# Basic run
python skills/openpiv/scripts/runner.py \
    --image img1.bmp --image img2.bmp --output_dir results --verbose

# Tuned parameters with dynamic masking
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	フレーム間の時間（s）
--scaling	96.52	スケーリング係数。物理単位あたりのピクセル数（例：px/mm）
--threshold	1.05	peak2peak の信号対雑音比のしきい値
--mask	none	none または dynamic（openpiv.preprocess.dynamic_masking）
--mask_method	intensity	edges または intensity。--mask dynamic と併用する場合のみ使用
--drop_invalid	オフ	補間値を残す代わりに、フラグが付いたベクトルを NaN にする
--verbose	オフ	進行状況のメッセージを表示

OpenPIV に同梱されている画像ペアを使って、インストールが一通り正常に動作することを確認してください。

python skills/openpiv/scripts/run_example.py --output_dir /tmp/openpiv-demo

出力ファイル

• vectors.txt — タブ区切り、%.4e 形式、# x y u v flags mask のコメントヘッダー付き
• params.npz — x、y、u、v、flags 配列を含む NumPy アーカイブ
• vector_field.png — 1 枚目のフレーム上に描画したベクトル場

# x	y	u	v	flags	mask
2.1757e-01	3.5226e+00	-6.2220e-02	-2.7081e+00	0.0000e+00	0.0000e+00
4.8695e-01	3.5226e+00	-3.1587e-01	-2.9800e+00	0.0000e+00	0.0000e+00

flags は浮動小数点数として書き込まれ、有効なベクトルは 0、フラグ付きのベクトルは 1 です。

ベストプラクティス

パラメーターの選択

1. ウィンドウサイズ — ほとんどの場合、32×32 が適しています。より粗い解像度で相関を高めるには 64/128、ノイズを許容してより細かい解像度にするには 16/24 を使用します。
2. オーバーラップ — ウィンドウサイズの 50～75％です。
3. しきい値 — より多くのベクトルを除外するには値を上げます。sig2noise_method を切り替えた後は、必ず再調整してください。
4. スケーリング係数 — キャリブレーショングリッドなどの既知の基準に照らして校正し、単位を正しく扱ってください（OpenPIV の test1 チュートリアルデータにおける 96.52 は px/mm です）。

画像品質

• 粒子が見やすく均等に分布しており、1 つの検査ウィンドウあたり 5～10 個ある
• 飽和または露出過多の領域がない
• 背景ノイズが最小限であること。実行全体を通じて背景差分を検討してください

処理のヒント

1. まずデフォルト設定から始め、得られたベクトル場を見ながら調整します。
2. s2n の分布を確認してください。中央値が低い場合は、しきい値ではなく相関が不良であることを示しています。
3. 早い段階で可視化してください（一様なベクトルやエッジアーティファクトなど、明らかな問題をすぐに確認できます）。
4. 速度勾配や変位が大きい流れには、複数パス（windef）を使用してください。
5. 反射や固体境界は、ベクトルを生成させるのではなく、マスクしてください。

リソース

references/

• advanced_algorithms.md — 相関およびサブピクセル法、複数パスのウィンドウ変形、PIVSettings のフィールド、3D および位相分離モジュール

アルゴリズムや設定の詳細情報が必要な場合は、リファレンスを読み込んでください。

関連スキル

• scvelo：スプライス層と非スプライス層（velocyto、ST…由来）を持つ単一細胞 RNA-seq AnnData オブジェクトに対して、scVelo を使用して RNA 速度解析を実行します。