이 가이드는 일본어만 있습니다(사람이 쓴 산문이라 기계적으로 치환하지 않았습니다).
2 枚の画像から、どこがどれだけ動いたかを測る層です。入力はトレーサ粒子を写した画像対、出力は窓ごとの変位ベクトル場 (2, h, w)(成分 (dy, dx)、単位は画素/フレーム)。流体計測(PIV)が本来の用途ですが、原理は相互相関なので、模様があるものなら粒子でなくても動きが取れます。
26 op / 7 カテゴリ(numpy と scipy のみ。台帳は opspiv.py、実体は pivops.py と dic.py):
• synth(3) — piv_synth_particles / piv_synth_pair / piv_synth_sequence: 既知の変位場を持つ画像対と画像列を作る。この族の全テストの真値の供給源。
• estimate(4) — piv_cross_correlate / piv_multipass / piv_deform_pass / piv_ensemble_correlate: 本体。窓ごとの FFT 相互相関、粗→細の多段、窓変形、相関の合算。
• validate(2) — piv_outlier_mask / piv_replace_outliers: 正規化中央値検定と穴埋め。
• field(8) — piv_vorticity / piv_divergence / piv_flow_magnitude / piv_to_velocity / piv_velocity_gradient / piv_q_criterion / piv_swirling_strength / piv_strain_rate: 場の微分・不変量・単位変換。
• visualise(2) — piv_flow_to_rgbimage / piv_line_integral_convolution: 出口。これが無いと flow2d は「作れるが見られない」型になる。
• solid(3) — strain_from_displacement / correlation_quality / speckle_quality: 固体側(DIC、デジタル画像相関)。同じ相関器を使うが、出す量が違う(下の節)。
• assess(4) — piv_sample_at_windows / piv_error_stats / piv_peak_locking / piv_time_statistics: 真値との突き合わせ、系統誤差、時間統計。
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 本だけでした。
Mermaid 図(ソース):
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 つを残しているのは選択肢のためではなく、系統誤差の違いを測れるようにするためです。
| 変位場 | 閉形式 | 渦度の実測 | 発散の実測 |
|---|---|---:|---:|
| 剛体回転 ω=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)) は画面上では時計回りに見える場)。規約を使って書いたテストでは規約の反転を捕まえられないので、向きが分かる最小の場で別に固定してあります。
flow2dflow_dense の述語は (3, D, H, W) 限定なので、2 成分の平面フローはそもそも該当しません。名前を借りると台帳が「3 成分を返す」と宣言しながら 2 成分を返すことになります。
型を増やすときの本 repo の条件(種を持つ op が無ければ永久に未実行になる)は満たしています —— 生成が 7 op、消費が 13 op、族の中で閉じています。出口(可視化 2 op)も必ず持たせています —— 作れるが見られない型は連鎖の途中で行き止まりになります。両方向の fail-closed も実測済み:
• reprconv.flow_magnitude(2-D flow) → ValueError: this op takes DENSE scene flow (3, D, H, W) …
• piv_vorticity(3-D scene flow) → ValueError: flow must be (2, h, w) …
型では守れないものも正直に書いておきます。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 引数 | 既定値なし | 既定値を置くと単位事故が既定になる |
• 窓の大きさ —— 粒子が窓あたり 5-10 個になるように選ぶ。密度 0.02 なら 32x32 で約 20 個。窓 16 まで下げると 5 個になり、実測で RMS が 0.060 → 0.113 と悪化した(多段 64→32→16)。細かくすれば良くなるわけではない。
• 多段 —— 格子を密にしたいときに効く。実測では窓 64 単段(格子 7x7、RMS 0.038)に対して 64→32(格子 15x15、RMS 0.060)。密度と精度は交換関係であって、多段は精度を上げる道具ではなく、精度をあまり落とさずに密度を上げる道具。
• 窓関数 —— hann は縁のリークを減らすが有効粒子も減らす。実測(補正なし・窓 32)で RMS 0.154(hann)対 0.213(none)なので既定は hann。
• **peak_ratio** —— 第 1 ピーク / 第 2 ピーク。1 に近い窓は当てにならない。雑音を足すと下がることを実測で固定してある。
渦度だけを見るとせん断層も光ります。層流の壁近傍が渦のように見えるのはそのためで、渦を取り出したいときは速度勾配テンソルの不変量を使います。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 だけでは薄い領域まで拾うので、閾値をどう決めたかを書かない渦可視化は、絵の美しさが閾値の産物である可能性を隠しています。閾値を振ったときの面積変化を併記してください。
• piv_flow_to_rgbimage — 色相 = 向き、明度 = 速さ。reprconv.flow_to_rgbimage の 2 次元版(あちらは 3-D シーンフロー専用)。**scale を省くと図ごとに色の意味が変わる**ので、複数の図を並べるなら固定すること。色相環の凡例は図の側で必ず一緒に焼く。
• piv_line_integral_convolution — 流れに沿って白色雑音を平均した模様。矢印は密にすると潰れ、疎にすると構造を見落としますが、LIC は画素ごとに積分するので密度の選択が要りません。向きの情報は落ちる(前後を区別しない)ので、回転の向きを見たいときは矢印か色相図と併用します。積分長を伸ばすと滑らかになりますが、渦の芯のように曲率が大きい場所では構造が伸びて嘘になります。
整数ずらしは「窓の中で変位が一定」を仮定します。回転やせん断ではこれが崩れて相関ピークが潰れるので、予測の半分ずつ 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 はコマごとに場全体を揺らす形にしてあります。
公開データ(SPLEEN 高速タービン翼列 PIV、Zenodo 10256759、CC BY 4.0、27 MB、MATLAB v5 の .mat なので scipy.io.loadmat がそのまま読む)で場の op を回した記録。
一致したもの:
• 流れ角 —— 同梱の beta と atan2(V_tan, V_axi) の差が中央値 0.0006 度(座標系の解釈が正しい確認)。
• 音速 —— V_mag_2D / Mach_2D = 327.2 m/s。
• 速さ —— 自前の piv_flow_magnitude と同梱値の差が中央値 0.067 m/s(約 296 m/s に対して相対 2e-4)。
ここで一度、失敗した報告をしました。 最初に出した数字は「最大差 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 を伝播させるので、埋めるかどうかは呼ぶ側の判断です。
流体の 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 を必須にしてあり、
• "infinitesimal" = ∂u/∂x(教科書の微小ひずみ)。回転で cosθ-1 を返す。
• "green" = Green-Lagrange ∂u/∂x + ½((∂u/∂x)² + (∂v/∂x)²)。**剛体回転で
代数的に厳密 0**(0.5〜5° で 5.3e-14 〜 6.5e-13)。
既定は置きません。実測では 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 のデスペックルで、これとは別物です。
• 公開データ —— 合成で閉じたら RSPID(合成粒子画像 + 真値)、円柱後流の時間分解 PIV(ストローハル数 0.2 前後という独立検算が効く)、タービン翼列。RAD コーパス fullseye_poc_datasets_corpus_v2 のノート 060-064 に URL の実測結果とライセンスの注意がある。
• 航空・流体 —— 渦度から Q 基準・λ2 へ(速度勾配テンソルの不変量)。同コーパスのノート 036。
• 既存の族 —— 出力の image2d は 2-D op 一式(閾値・morphology・疑似カラー・図注)へそのまま渡せる。時系列にすれば videostream 族。