openpiv

OpenPIV

概要

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

対象は openpiv 0.26.1 です。合成データによる検証は Python 3.13、NumPy 2.5.3、SciPy 1.18.1、scikit-image 0.26.0、Matplotlib 3.11.2 で実行しました。Rust はインストールされていません。検証済みの経路を使用するには backend="scipy" を指定してください。詳細はレビューの根拠を参照してください。

使用する場面

このスキルは、実験的な 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.26.1"

画像ペアに対して 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")

# Always returns (u, v, s2n); s2n is all NaN if sig2noise_method=None.
u, v, s2n = pyprocess.extended_search_area_piv(
    frame_a.astype(np.float32),
    frame_b.astype(np.float32),
    window_size=32,
    overlap=12,
    dt=0.02,
    search_area_size=38,
    correlation_method="linear",   # circular also supports extended search
    normalized_correlation=True,
    backend="scipy",
    sig2noise_method="peak2peak",
)

x, y = pyprocess.get_coordinates(
    image_size=frame_a.shape,
    search_area_size=38,
    overlap=12,
    center_on_field=False,         # matches sliding-window positions
)

# flags is a boolean array: True marks a spurious vector.
flags = validation.sig2noise_val(s2n, threshold=1.05) | ~np.isfinite(s2n)
flags |= ~np.isfinite(u) | ~np.isfinite(v)
u, v = filters.replace_outliers(u, v, flags, method="localmean", max_iter=3, kernel_size=2)

# Scale, then convert to right-handed image-boundary coordinates.
x, y, u, v = scaling.uniform(x, y, u, v, scaling_factor=96.52)
y, v = frame_a.shape[0] / 96.52 - y, -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%）。overlap を大きくするとベクトル密度と計算コストが増加しますが、隣接するベクトルは独立ではなく相関を持つようになります。

search_area_size — 2 枚目のフレームで探索する窓。window_size 以上でなければなりません。数ピクセル大きくすると、より大きな変位に対応できます。"circular" と "linear" の両方が拡張探索をサポートします。CLI はこれに対して normalized_correlation=True を指定したゼロパディングの "linear" を選択しますが、これは任意の変位を復元したり、偽のピークを排除したりするものではありません。グリッドのストライドは search_area_size - overlap なので、overlap は探索領域より小さくなければなりません。

経験則: 最大変位は window_size のおおよそ 4 分の 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.preprocess にあり、(image, mask) を返します。符号なし減算によるアーティファクトを避けるため、浮動小数点のコピーを使用してください。解析前に両方のフレームでマスクを確認してください。

from openpiv import preprocess

# Intensity masking uses Otsu; the threshold argument is for edges.
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_a | mask_b を除外する物理領域として保持してください。CLI はこれを検査中心でサンプリングし、保存し、置換後もそれらの速度を NaN のまま保持します。オブジェクトにまたがる窓は依然としてバイアスを受ける可能性があります。実験に応じてマスクを検査・膨張してください。0.26.1 では method="edges" は bool ではなく uint8 でインデックスするため、画像が破損したり失敗したりする可能性があります。CLI はこれを拒否します。強度マスキングを使用するか、カスタムワークフローで検証済みの boolean マスクを使用してください。

マルチパス処理

マルチパス（ウィンドウ変形）は openpiv.windef にあり、PIVSettings データクラスによって駆動されます。 pyprocess にはマルチパスのエントリポイントはありません。

import numpy as np
from openpiv import scaling, windef

settings = windef.PIVSettings(
    windowsizes=(64, 32, 16), overlap=(32, 16, 8), num_iterations=3,
    backend="scipy", sig2noise_method="peak2peak", sig2noise_threshold=1.05,
)

x, y, u, v, flags = windef.simple_multipass(
    frame_a.astype(np.float32), frame_b.astype(np.float32), 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 をゼロで埋め、座標を変換します。 そのフラグには注意してください。修復済み・ゼロ埋め済みの出力は独立した計測値ではありません。 必ず settings を渡してください。settings を渡さない経路は 0.26.1 ではウィンドウサイズが 2 つである一方、反復回数が 3 回になっています。

単位の罠: simple_multipass とそのエイリアス multigrid_windef は settings.dt と settings.scaling_factor を無視します。これらの配列は上記のように変換してください。バッチの windef.piv(settings) は保存前に両方を適用します — そのファイルを再度スケーリングしないでください。simple ラッパーはすべてのバッチ前処理・出力スイッチを尊重するわけではありません。

個々のパスを制御するには、windef.first_pass と windef.multipass_img_deform が より低レベルな構成要素です。

検証と後処理

検証メソッド

すべてのバリデータは、True が疑わしいベクトルを示すブール配列を返します。また、非有限の u、v、s2n も 拒否してください。sig2noise_val だけでは NaN を検出できません。

# 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/フレームの変位は 150 px/s として渡されます。上記のしきい値はそのケースに適しています。PIV の文献や PIVSettings.min_max_u_disp が用いる (-30, 30) という値は px/フレームの制限であり、これを px/s 出力に適用すると フィールド全体が拒否されてしまいます。変換前に変位を検証するには dt=1 で相関を取るか、実際の速度単位（時間とキャリブレーション係数の両方を含む）にしきい値を調整してください。

外れ値の置換

u, v = filters.replace_outliers(
    u, v, flags, method="localmean", max_iter=3, tol=1e-3, kernel_size=2
)

method は "localmean"、"disk"、"distance" を受け付けます。無効な値を指定すると、欠損セルを埋める際に ValueError が発生します。"distance" には kernel_size>=2 を使用してください。0.26.1 ではサイズ 1 はすべての近傍重みをゼロに切り詰め、穴を残します。置換はフラグが付いた 位置を補間値で埋めます — その後でそれらを 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/Inf は欠測観測値です。入力が変更される可能性があるため、コピーを渡してください。先に穴をゼロ埋めすると、それらを計測されたゼロとして誤って扱ってしまいます。

from openpiv.smoothn import smoothn

u_smooth, *_ = smoothn(np.asarray(u, dtype=float).copy(), s=0.5)  # larger == smoother
v_smooth, *_ = smoothn(np.asarray(v, dtype=float).copy(), 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=False,
)
image = tools.imread("frame_a.bmp")
height, width = image.shape
ax.imshow(image, cmap="gray", origin="upper", zorder=-1,
          extent=(0, width / 96.52, 0, height / 96.52))
fig.savefig("vector_field.png", dpi=150, bbox_inches="tight")
plt.close(fig)

画像の範囲を明示的に設定してください。このヘルパーの on_img=True は最後のベクトルから範囲を推測し、 ウィンドウが未使用のマージンを残す場合にオーバーレイを引き伸ばすことがあります。CLI は実際の画像境界でプロットします。

カスタム可視化

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()             # excludes flagged and masked vectors
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)

シングルパスの拡張サーチグリッドでは、間隔は物理単位で (search_area_size - overlap) / scaling_factor となります。これが (window_size - overlap) / scaling_factor に簡約されるのは、2 つのウィンドウサイズが一致する場合のみです。 保存された x 座標と y 座標の差分を使用することを推奨します。特にマルチパス処理後は重要です。 dx=1.0 のままにすると、単位長さあたりではなくグリッドセルあたりの渦度となります。 OpenPIV の座標生成を参照してください。

符号の規約: CLI は y = image_height / scaling_factor - y を書き込み、v を符号反転し、 行は画像の順序のまま保持します。保存された y は配列の下方向に向かって減少します。単体の形式では 符号付きの dy(ここでは負)が必要です。そうでなければ PIVAnalyzer を使用してください。これは向きを読み取ります。正の dy は y 微分の寄与を変え、実際の回転を打ち消してしまうことさえあります。

ひずみ速度

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 フレームの空間平均を差し引くと空間分散が測定されますが、これは 均質性・エルゴード性の仮定が正当化される場合にのみ乱流強度と等しくなります。レイノルズ分解には アンサンブルが必要です。各実現からその平均場を差し引いてください。レガシーの tke キーは 2 成分の空間分散エネルギーのみであり、完全な乱流運動エネルギーではありません。アナライザーはデフォルトで フラグ付きベクトルを除外します。compute_statistics(include_interpolated=True) は修復されたものを含めます。

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	0.26.1 では intensity のみ使用可能です。edges は拒否されます
--backend	scipy	scipy、auto、または rust。Rust は拡張機能が必要です
--drop_invalid	off	補間された値を保持する代わりに、フラグ付きベクトルを NaN にします
--verbose	off	進捗メッセージを表示します

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、mask、s2n、タイミング/キャリブレーション/ウィンドウのパラメータ、 要求されたバックエンド、および OpenPIV のバージョン。座標は物理単位を使用し、速度は単位/秒を使用します。
• vector_field.png — 最初のフレーム上に描画されたベクトル場

# 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 は float で、0 は有効、1 はフラグ付きを表します。mask=1 は除外されたオブジェクト領域を フラグとは独立に示します。測定されたフィールドとともに両方を保持してください。

ベストプラクティス

パラメータの選択

1. ウィンドウサイズ — 32×32 がほとんどのケースに適しています。64/128 はより粗い解像度で より良い相関を得る場合に、16/24 はノイズと引き換えに、より細かい解像度を得る場合に使用します。
2. 重なり — ウィンドウサイズの 50～75%。
3. しきい値 — より多くのベクトルを棄却するには上げます。sig2noise_method を切り替えた後は 必ず再調整してください。
4. スケーリング係数 — キャリブレーショングリッドなどの既知の基準に対してキャリブレーションを行い、 単位を揃えてください(OpenPIV の test1 チュートリアルデータの 96.52 は px/mm です)。

画像品質

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

処理のヒント

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

リソース

references/

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

詳細なアルゴリズムや設定情報が必要な場合は、このリファレンスを参照してください。