fullseye

粒子画像流速測定(画像対から密な変位を測る) — 使い方ガイド

この族は何をする道具箱か

2 枚の画像から、どこがどれだけ動いたかを測る層です。入力はトレーサ粒子を写した画像対、出力は窓ごとの変位ベクトル場 (2, h, w)(成分 (dy, dx)、単位は画素/フレーム)。流体計測(PIV)が本来の用途ですが、原理は相互相関なので、模様があるものなら粒子でなくても動きが取れます。

26 op / 7 カテゴリ(numpy と scipy のみ。台帳は opspiv.py、実体は pivops.py と dic.py):

なぜ足したか —— 在庫を数えた結果

2026-09-06 に型付きカタログ全体(900 op)を走査したところ、画像対から密な変位を出す op が 1 つも無いことが分かりました。近いものはあります:

既にあるもの 何をするか なぜ代わりにならないか
scene_flow_lk Lucas-Kanade の密なフロー 3 次元の体積用((3,D,H,W))。平面の画像対は入らない
estimate_flow / nearest_neighbor_flow シーンフロー 点群同士。画像を食わない
correlation_score 相関 体積同士のスカラ 1 個。場を返さない
match_keypoints 対応点 疎。窓ごとの密な場ではない
optical_flow_magnitude_stream 動画の流れの大きさ 大きさのみ。ベクトルを返さない

flow_dense 型は存在しますが、述語が ndim == 4 and shape[0] == 3 で、作る op は scene_flow_lk の 1 本だけでした。

使う順序

flowchart TD
    S["piv_synth_pair<br/>既知の変位場 → 画像対 + 真値"]
    IMG["粒子画像 a, b<br/>(H, W)"]
    S --> IMG
    REAL["実験・公開データ<br/>(RSPID, 円柱後流 …)"] --> IMG
    IMG --> CC["piv_cross_correlate<br/>窓ごとの FFT 相互相関"]
    IMG --> MP["piv_multipass<br/>粗→細(予測変位つき)"]
    CC --> F["flow2d (2, h, w)<br/>(dy, dx) [px/frame]"]
    MP --> F
    F --> OM["piv_outlier_mask<br/>正規化中央値検定"]
    OM --> RO["piv_replace_outliers"]
    RO --> F2["整えた flow2d"]
    F2 --> VO["piv_vorticity<br/>渦度"]
    F2 --> DI["piv_divergence<br/>発散(非圧縮なら 0 = 独立検算)"]
    F2 --> MG["piv_flow_magnitude"]
    F2 --> VE["piv_to_velocity<br/>m/s(画素寸法と Δt が必須)"]
    S --> TR["真値 (2, H, W)"]
    TR --> SA["piv_sample_at_windows<br/>窓格子へ落とす"]
    SA --> ES["piv_error_stats<br/>偏りと散らばりを分ける"]
    F2 --> ES
    F2 --> PL["piv_peak_locking<br/>系統誤差の強さ"]

最小の例(そのまま動きます)

import numpy as np

import pivops

# 1) 既知の変位場を決めて画像対を作る —— 真値は定義そのもの
cy, cx = 127.5, 127.5


def rigid_rotation(rows, cols):
    """画面上で時計回りに回る場。渦度は -2w になる(下の assert 参照)。"""
    w = 0.01
    return w * (cols - cx), -w * (rows - cy)


a, b, truth = pivops.piv_synth_pair((256, 256), rigid_rotation,
                                    density=0.02, diameter_px=2.5, seed=7)

# 2) 粗い窓から細かい窓へ(多段)
flow, info = pivops.piv_multipass(a, b, windows=(64, 32), overlap=0.5)

# 3) 外れ値を見つけて埋める
bad = pivops.piv_outlier_mask(flow, threshold=2.0)
flow = pivops.piv_replace_outliers(flow, bad, "median")

# 4) 場の量
vort = pivops.piv_vorticity(flow, info["step"])
div = pivops.piv_divergence(flow, info["step"])

# 5) 真値と突き合わせる(格子を揃えてから)
t = pivops.piv_sample_at_windows(truth, info)
stats = pivops.piv_error_stats(flow, t)

assert stats["rms"] < 0.12                       # 実測 0.05 前後
assert abs(stats["bias_dy"]) < 0.02              # 零方向への偏りは補正済み
# 渦度は閉形式 -2w = -0.02 に一致する(**負**なのが規約どおり)
assert abs(np.mean(vort[1:-1, 1:-1]) - (-0.02)) < 0.0015
# 発散は 0 —— 真値と比べるのとは**別経路**の検算
assert abs(np.mean(div[1:-1, 1:-1])) < 5e-4

# 6) 物理速度へ。画素寸法と時間差は**必須引数**(既定値を置いていない)
v = pivops.piv_to_velocity(flow, pixel_size_m=1e-5, dt_s=2e-4)   # 10 um/px, 200 us
assert np.isfinite(v).all()

零方向への偏りと、その補正(この族でいちばん大事な実測)

素の相互相関は変位を零へ引き寄せます。窓をずらすと重なる領域が減り、相関の値そのものが変位とともに落ちるからです。実測(窓 32、Hann、dx を振る):

真の dx [px] 偏り(補正なし) 偏り / (d/N)
0.5 -0.0203 1.30
1.0 -0.0402 1.29
2.0 -0.0802 1.28
4.0 -0.1600 1.28
6.0 -0.2406 1.28

比が一定であることが、原因の説明(重なり面積)が合っている証拠です。値が小さいことではありません。

補正は 2 つを必ず対で入れます: (a) 窓関数の自己相関で割る(normalize="overlap")、(b) 探索を窓の 1/4 に絞る(search_limit=0.25、PIV の「1/4 則」)。片方だけだと悪化します:

設定 偏り dy 偏り dx RMS
両方あり(既定) -0.0045 +0.0039 0.0648
正規化なし -0.0867 +0.1030 0.1535
正規化だけ(探索無制限) -0.1910 +0.2329 2.9094

正規化だけを入れると RMS が 45 倍に悪化します。縁では割る量が 0 に近づき、そこに偽のピークが立つためです。「補正を足したのだから良くなったはず」を測らずに信じない、という戒めがそのまま設定の既定になっています。

既知の系統誤差 —— 出るはずのものが出るか

サブピクセル推定は、真の変位の小数部を整数へ引き寄せる偏り(ピークロッキング)を持ちます。小数部を 0 から 0.9 まで振った実測:

推定法 小数部誤差の RMS 最大絶対誤差 0.1 の答え 0.9 の答え
gauss3(既定) 0.0037 px 0.0088 px 0.099 0.896
parabolic 0.0104 px 0.0151 px 0.091 0.904
centroid 0.2259 px 0.3701 px 0.021 0.978

centroid の S 字は教科書どおりで、出ないほうがおかしい。3 つを残しているのは選択肢のためではなく、系統誤差の違いを測れるようにするためです。

閉形式との一致(320x320、多段 64→32)

変位場 閉形式 渦度の実測 発散の実測
剛体回転 ω=0.01 渦度 -2ω = -0.02、発散 0 -0.01997 -0.00006
一様膨張 s=0.01 発散 2s = 0.02、渦度 0 +0.00009 +0.01985
単純せん断 g=0.02 渦度 g = 0.02、発散 0 +0.01996 -0.00003

回転の渦度が負なのは規約どおりです。テストを書いたとき最初に +2ω と書いて落ちましたが、間違っていたのはテストの側でした((ω(c-cx), -ω(r-cy)) は画面上では時計回りに見える場)。規約を使って書いたテストでは規約の反転を捕まえられないので、向きが分かる最小の場で別に固定してあります。

型を 1 つだけ新設した理由 —— flow2d

flow_dense の述語は (3, D, H, W) 限定なので、2 成分の平面フローはそもそも該当しません。名前を借りると台帳が「3 成分を返す」と宣言しながら 2 成分を返すことになります。

型を増やすときの本 repo の条件(種を持つ op が無ければ永久に未実行になる)は満たしています —— 生成が 7 op、消費が 13 op、族の中で閉じています。出口(可視化 2 op)も必ず持たせています —— 作れるが見られない型は連鎖の途中で行き止まりになります。両方向の fail-closed も実測済み:

型では守れないものも正直に書いておきます。piv_to_velocity の出力は単位が m/s に変わりますが型は同じ flow2d です。単位を型で分ける案は、m/s を作る op が 1 本・消費する op が 0 本の型を生むので採りませんでした(「証拠が出てから増やす」の順序に反する)。代わりに画素寸法と時間差を必須引数にしてあります。

規約(取り違えると静かに間違う)

規約 値 破ったときに何が起きるか
変位の単位 画素/フレーム m/s と混ぜると桁が変わる。例外は出ない
成分順 (dy, dx) 転置した場が出る。絵は自然に見える
dy の正の向き 行が増える向き(画像の下) 上下反転
渦度の定義 d(dx)/dy - d(dy)/dx 反時計回りが正。符号を決めずに書くと渦の向きが黙って反転する
ベクトルの位置 窓の中心(info["rows"] / ["cols"]) 真値と半窓ずれて比較され、勾配のある場で誤差が水増しされる
piv_to_velocity の 2 引数 既定値なし 既定値を置くと単位事故が既定になる

使い分けの目安(実測に基づく)

渦とせん断を分ける(派生 op)

渦度だけを見るとせん断層も光ります。層流の壁近傍が渦のように見えるのはそのためで、渦を取り出したいときは速度勾配テンソルの不変量を使います。3 つとも同じ 4 つの微分から出るので、piv_velocity_gradient が一度にまとめて返します。

解析場を直接置いた実測(PIV を通さない、定義そのものの検算):

場 渦度 Q 基準 渦回転強度 λ_ci ひずみ速度
剛体回転 ω=0.01 -0.020000 +1.0e-4 (= ω²) 0.010000 (= ω) 0.000000
一様膨張 s=0.01 0.000000 -1.0e-4 (= -s²) 0.000000 0.020000 (= 2s)
単純せん断 g=0.02 +0.020000 0.000000 0.000000 0.020000 (= g)

最下行が要点です。せん断の渦度は回転と同じ大きさなのに、Q も λ_ci も 0。この差がこの 2 つを足した理由そのものです。

Q も λ_ci も閾値を必要とする量です。Q > 0 だけでは薄い領域まで拾うので、閾値をどう決めたかを書かない渦可視化は、絵の美しさが閾値の産物である可能性を隠しています。閾値を振ったときの面積変化を併記してください。

見せ方(可視化の 2 つ)

精度を上げる 2 つの道(実測つき)

窓変形 —— 窓の中で変位が変わる場に効く

整数ずらしは「窓の中で変位が一定」を仮定します。回転やせん断ではこれが崩れて相関ピークが潰れるので、予測の半分ずつ 2 枚を逆向きに歪めてから相関します(中央差分の変形)。

手順 RMS 偏り (dy, dx)
多段 64→32 のみ 0.0576 (-0.0055, -0.0014)
+ 窓変形 1 段 0.0272 (+0.0012, -0.0003)

回転場(ω=0.01、256x256)で 2.1 倍。

アンサンブル相関 —— 疎で雑音の多いデータに効く

1 対ずつ測って平均すると「外れたベクトルの平均」になります。相関マップの段階で足すと、弱いピークが同じ場所に積み上がって立ち上がります(定常流が前提)。

手順 RMS 1 px 超の外れ
1 対のみ 4.21 47 / 81
10 対を合算 1.55 9 / 81

密度 0.004(窓 32 あたり約 4 個)+ 雑音 0.25 という厳しい条件での実測です。完全には直りません —— 2.7 倍良くなるが、9 本はまだ外れたまま。そう書いておきます。

時間統計 —— 平均・変動・レイノルズ応力

piv_time_statistics は画像列から、時間平均・変動の RMS・レイノルズ応力 ⟨u’v’⟩ を返します。乱流の記述はこの 3 つが出発点で、⟨u’v’⟩ の符号と大きさが運動量輸送そのものになります。

検証は piv_synth_sequence に既知の乱れを仕込んで行います。実測: 仕込み 0.3 px → 測定 0.287 px、定常流(乱れ 0)なら 0.013 px。

★ ここで一度失敗しています。最初は独立な画像対を並べたものを「列」として渡し、平均 0.73 px に対して変動の RMS が 4.05 px という数字を得ました。隣り合う 2 枚に対応関係が無いので、流れではなく列の作り方を測っていたのです。piv_synth_sequence は同じ粒子を追い続けます。

もう一つ。乱れを粒子ごとに独立に振ると、窓の中で平均されて粒子数の平方根ぶん小さくなり、仕込み 0.3 に対して 0.15 しか出ませんでした。これは PIV の性質であって誤りではありませんが、時間統計の検証には使えないので、jitter はコマごとに場全体を揺らす形にしてあります。

実データに当てたときに起きたこと(2026-09-06)

公開データ(SPLEEN 高速タービン翼列 PIV、Zenodo 10256759、CC BY 4.0、27 MB、MATLAB v5 の .mat なので scipy.io.loadmat がそのまま読む)で場の op を回した記録。

一致したもの:

ここで一度、失敗した報告をしました。 最初に出した数字は「最大差 200 m/s、相対 0.81」でした。事実ですが、分布はこうです:

差の絶対値 50 % 90 % 99 % 最大
[m/s] 0.067 1.099 15.864 200.144

中央値と最大で 3 桁半違います。最大値を見出しにすると「合っていない」に、中央値を見出しにすると「合っている」に見える —— どちらも本当なので、両方出すのが正しい報告です。

食い違いの原因は実装ではなく比べている量の違いでした。V_mag_2D は |平均速度| ではなく 平均した速さ で、イェンセンの不等式から前者以上になります。裏づけ: 差と (STD_Vaxi² + STD_Vtan²) / (2|V|) の相関 r = 0.879、差が 20 m/s を超える点の乱流強度は中央値 42 %(全体 2.2 %)—— つまり後流の中に集中していました。

公開された場と突き合わせるときは、量の定義を先に一致させること。名前(V_mag_2D)は平均の取り方を書いていません。

もう一つ実データで効いた事実: 欠測が 67.5 % ありました(翼の影と視野外)。nan を 0 で埋めて微分すると、翼の輪郭に沿って偽の勾配が立ちます。この族の場の op は nan を伝播させるので、埋めるかどうかは呼ぶ側の判断です。

固体を測る(DIC)—— 同じ相関器、違う出口

流体の PIV と固体の DIC は同じ窓相関を使います。違うのは出す量で、流体は 速度と渦度、固体はひずみ。ここを混ぜると静かに間違うので、2026-09-06 に examples/poc_dic_strain.py で実測してから 3 op だけ足しました。

piv_strain_rate を固体に使ってはいけない

docstring には「剛体回転では 0 になる」と書いてあり、流体としては正しい —— 線形化した回転 u = -ωy, v = ωx では確かに 0 です。ところが DIC が測るのは 有限回転で、そちらでは

θ piv_strain_rate 真のひずみ
0.5° +76 µε 0
1.0° +305 µε 0
2.0° +1218 µε 0

鋼の降伏ひずみが約 2000 µε なので、2 度傾いただけで降伏の 6 割に相当する 嘘のひずみが出ます。strain_from_displacement(u, v, window, method) は method を必須にしてあり、

既定は置きません。実測では Green の補正項が勾配の推定誤差を拾うので、 勾配の散らばりが大きい推定器では補正が過剰になって悪化することがあります (PoC の 6 節に表があります)。どちらを使うかは呼び手が決めるべき分岐です。

微分は窓最小二乗で

piv_velocity_gradient は np.gradient(2 点差分)です。同じ流れ場・500 µε で

方法 平均 散らばり
piv_velocity_gradient 500.0 µε 136.9 µε
strain_from_displacement(窓 9) 500.0 µε 12.4 µε

平均は同じで散らばりが 11 倍違います。領域平均で使うなら差は出ませんが、 1 点の値を見るなら効きます。

品質で切る —— peak_ratio では足りない

piv_cross_correlate の info["peak_ratio"] は相関ピークの鋭さですが、 対応が無くなった領域を分離できません。60×60 の領域を別の絵に貼り替えた 実測で、内 1.184 / 外 1.329 —— 分布が重なります。ZNCC 係数なら 0.100 / 0.999。

切り方 変位の RMS
切らない 1.9122 px
peak_ratio >= 1.2 1.0952 px(1.7 倍)
correlation_quality >= 0.8 0.0036 px(537 倍、82.8 % を残す)

correlation_quality(ref, cur, flow, info) は相関器ではありません —— 既にある変位場を採点するだけなので、piv_cross_correlate でも optical_flow_lk でも demons_register の結果でも使えます。

撮る前にスペックルを測る

speckle_quality(img) は平均輝度勾配(MIG)・被覆率・平均斑点径を返します。 PoC の 9 節のとおり、同じ推定器でも斑点の粒径で偏りが 20 倍変わります (1σ 0.6 px で 0.0113 px、1σ 4.0 px で 0.0005 px)。撮ってから気づくと撮り直しです。 speckle_filter は SAR のデスペックルで、これとは別物です。

次にどこへ繋がるか

This site is open source. Improve this page.