動いているのに見えない層です。カメラの前で機械の枠が共振で 0.1 画素だけ呼吸している、配管の壁が脈打っている、ボルト継手が 1 ミクロンずつ緩んでいる — どれも記録には入っているのに、人間の目にも普通の運動解析にも掛かりません。そこから欲しいものは 2 つあって、別の問題です。
どちらも同じ量から出ます。向きつき帯域成分の局所位相です。帯域制限された画像成分が d だけ平行移動すると、その成分の局所位相は -k·d だけずれます(k は成分の局所空間周波数、単位 rad/px)。位相はサブピクセル変位の線形な符号化なので、位相を α 倍すれば変位が α 倍になります — 運動場を一度も推定せずに、厳密に。9 op / 5 カテゴリ(numpy + scipy のみ、台帳は opsmotionmag.py、実体は motionmag.py):
synthesize_translation: 既知振幅・既知周波数のサブピクセル平行移動クリップ。フーリエ位相ランプで動かすので補間誤差がゼロで、真値が閉形式。この族の全主張はこの合成に対して測ってあります。complex_steerable_decompose / complex_steerable_reconstruct: 1 枚を向きつきの解析信号(複素サブ帯域)に割り、戻す。tight frame として組んであるので往復が機械精度で厳密です。ここが厳密でないと「増幅のせいか再合成のせいか」が永久に切り分けられません。temporal_bandpass / temporal_band_power / band_snr: 画素ごとの時間軸を帯域で切る、その帯域のパワー地図を返す(= どこが何 Hz で振動しているか)、そしてクリップの帯域内容と雑音床を分解して測る。motion_magnify: 帯域内の位相を alpha - 1 倍して足し戻し、再合成する。返り値は動画と SNR 一式。phase_displacement / displacement_series: 画素ごとの変位場 dx / dy(と観測可能性の rank)、および全画面の変位波形 (T, 2)。後者はそのまま dsp.spectrum に流せば共振周波数が読めます。alpha は変位利得として定義してあります: 出力の変位が alpha * d。alpha = 1 が恒等、0 が運動の除去、-1 が反転です。公開文献は増幅後を 1 + alpha_paper 倍と書きますが、本族の alpha は 1 + alpha_paper に当たります。理由は検証可能性で、呼び出し側が書いた数と測定値を、足し引きなしにそのまま比較できるようにするためです。
データ種は既存語彙の再利用が基本です: image2d(帯域パワー地図)、table(分解結果 dict・SNR 計測 dict・増幅結果 dict・変位場 dict)、pairs((n,2) 配列 = 変位波形。MTF 曲線や funct1d と同じ規約)。新語は video の 1 つだけで、(T, H, W) の float64 フレーム列です。これは videops が使っている規約をそのまま踏襲しています。
「どこがどれだけ動いたか」と「サブピクセルの周期運動を帯域で選んで増幅する」は別の問題です。 前者はオプティカルフローの守備範囲で、後者がこの族です。混同すると、動かない答えを何時間も追うことになります。
| やりたいこと | 使う op | 置き場所 |
|---|---|---|
| 1 画素以上の運動、独立に動く物体、遮蔽をまたぐ追跡 | optical_flow_lk / optical_flow_hs / warp_by_flow / track_points |
flow(輝度不変を空間窓で解く。この族は一切呼ばないし再実装もしない — レジームが違う) |
| 推定済みフロー場の解釈(全体運動モデル・残差・運動分割) | dominant_motion / residual_motion / motion_segments |
motion(無関係・不変) |
| 時間方向の低域平滑、背景モデル、運動エネルギー地図 | moving_average / spatiotemporal_gaussian / background_subtraction / motion_energy |
videops(temporal_bandpass は同じ族の帯域選択メンバであって、これらの重複ではない) |
| 方向づきエッジ応答(実数) | tf_steerable_filter |
backends_transform2(直交対でなく位相を持たず可逆でもないので増幅には使えない。別物としてそのまま残す) |
| FFT・複素画像・2-D 位相アンラップ・Wiener 復元 | cx_fft 系 / phase_unwrap / cx_wiener_deconvolve |
complexops(FFT 畳み込み・相関は filters_freq) |
| 変位波形の周波数解析 | spectrum / signal_features |
dsp(displacement_series の (T, 2) をそのまま食える。ラップし直さない) |
ステアラブル束だけは numpy.fft の上に自前で組んであります。tight frame でないと再合成が厳密にならないためで、厳密さがこのモジュールの契約そのものだからです。
全 op が入力を検証してから計算します。以下は 2026-09-01 の敵対監査で実際に見つかったバグか、それを塞ぐために書いた罠です。この族が探したのは「例外が出る」ことではなく「黙って間違った数字を返す」ことでした。
alpha = 2 が 2d ではなく 0.939 * (2d) を返し、alpha = 0 は運動の 7.9 % を消し残しました。恐ろしいのは誤差が d に依存しないことです — 微小変位で消えないので「線形性の限界」では説明がつかず、非線形の言い訳が通りません。もっともらしく間違った数字の教科書例で、対策は synthesize_translation の既定波長を (8.0, 16.0) とオクターブを分けることにして、狭帯域条件を合成の側で可視化しました。np.unwrap が雑音帯を乱歩に変える。 大きな運動に届かせようと位相を時間方向にアンラップすると、信号を持たないサブ帯域ではフレーム間位相が -pi 〜 pi の一様乱数になり、アンラップはランダムウォークになります。実測: 0.05 px の運動に sigma = 0.02 の雑音を載せたクリップで、増幅器が実際に適用した最大位相増分は 12.27 rad — 意図は 0.039 rad でした。例外も NaN も出ず、単に違う動画が出ます。対策はアンラップの撤去で、angle(z * conj(z_mean)) は構成上 pi で有界なので、雑音帯の寄与は利得前で最大 pi に抑えられます。届く範囲は代わりに正直に開示しました(後述)。band_snr を増幅後クリップにそのまま当てると、起きていない改善を報告する。 band_snr は帯域内の雑音床を帯域外のビンから推定しますが、増幅は帯域外に触りません。増幅後の動画に当てると「帯域内パワーだけが alpha^2 倍、雑音推定は据え置き」となり、実測 alpha = 2 で +6.86 dB の運動 SNR 改善を報告しました。録画に無かった確からしさが計算で生えることはないので、これは嘘です。「評価指標を出力に当てると嘘になる」典型で、対策として motion_magnify は利得を知っているので入力側の雑音床を基準に補正した motion_snr_out_db を返します。素の result["snr_out"]["motion_snr_db"] は使わないでください。linear_regime が間違ったものを測っていた。 適用位相増分の最大値で判定していたため、コントラストを持たない雑音帯に引きずられて常に False になっていました。現在はコントラスト加重 RMS(phase_shift_rms_rad)で判定し、最大値は「どこかで起きた最悪」として別に返します。RMS が利得に厳密比例することはテストで固定してあります。rank を画素ごとに返します。float("30") は成功してしまうので、fps / f_lo / f_hi / alpha の文字列は ValueError(未パースの設定値が周波数として通り抜けるのを止める)。bool も拒否(True == 1 の暗黙昇格は、fps なら 1 Hz のタイムベース、利得なら恒等という別物になる)。complex も拒否(虚部の無言切り捨て)。配列側は complex / masked / NaN / Inf をすべて ValueError。f_hi > fps/2 はそのクリップに存在しない時間周波数なので、黙って低い周波数に畳まずに拒否します。fps と帯域の取り違え(fps=4 に 3-5 Hz を要求)もこの検査で捕まります。synthesize_translation も Nyquist を超える運動周波数を拒否します。fps <= 0、f_lo >= f_hi、そして DC に触れる通過帯域(f_lo <= 0)。最後のものは静止位相 = 「シーンがどこに在るか」を増幅してしまい、巨大で完全に架空の変位を作ります。fps/T より狭いと DFT ビンを 1 本も含まず、フィルタは黙ってゼロを返します。ビン数 0 は ValueError で、メッセージに分解能を書きます。band_snr は帯域外ビンが 1 本も残らない場合も拒否します(雑音床を推定する材料が無いため)。T = 1、2x2 フレーム、ぼろぼろのフレーム列。どれも数字を作らず、静止なら入力をそのまま返します。+-100 dB の窓にクランプし、クランプしたこと自体を snr_clamped で返します。「float64 が溢れた」と「答えが無限大」は別の主張だからです。MAX_FRAMES(4096)、MAX_FRAME_PIXELS(2²²)、MAX_VIDEO_ELEMENTS(2²⁴、軽い時間 op 用)、MAX_PYRAMID_ELEMENTS(2²²、複素コピーを同時に数本持つ op 用 ≈ 400 MB)、MAX_SCALES(8)、MAX_ORIENTATIONS(16)、MAX_FILTER_ELEMENTS(2²⁴)、MAX_ALPHA(200)。上限が無ければ「300 フレームの 1024x1024」という無害に見える指定が 30 GB を要求します。SNR を同時に報告しない増幅は嘘です。 増幅は帯域内の位相を alpha 倍しますが、それは帯域内の運動と帯域内の雑音を同じ係数で倍にします。したがって:
alpha^2 で増える一方、静止シーンは増えないからです。実測(64x64 / 64 frame / 32 fps、0.2 px の 4 Hz 運動、sigma = 0.01 のセンサ雑音、帯域 3-5 Hz):
| alpha | image SNR [dB] | image 変化 [dB] | motion_snr_out [dB] | motion 変化 [dB] | band_power_ratio |
|---|---|---|---|---|---|
| 1 | 29.2574 | -0.0000 | 11.9404 | +0.0000 | 1.000000 |
| 2 | 24.4304 | -4.8270 | 11.6285 | -0.3119 | 0.934861 |
| 4 | 18.9039 | -10.3535 | 11.2270 | -0.7134 | 0.857626 |
| 8 | 13.7428 | -15.5146 | 9.7565 | -2.1839 | 0.628551 |
倍ごとにおよそ 5 dB の画像 SNR を払っています。増幅した帯域が雑音収支を支配しきると、代数の漸近値である 20*log10(2) = 6.02 dB / 倍に近づきます。単調減少は sigma = 0.002 / 0.005 / 0.01 / 0.02 / 0.05 のすべてで実測・テスト固定済みです(雑音が大きいほど最初の一段の損失は小さくなります。帯域がもともと雑音床より下だと、倍にしても全変動のごく一部しか倍にならないためで、これは欠陥ではありません)。
band_power_ratio は測定値 band_power_out / (alpha^2 * band_power_in) です。1.0 なら増幅は線形に効いており、不足分は位相変調が高調波に投げ捨てたエネルギーです。仮定ではなく毎回測っているので、線形領域を出たことが返り値から分かります。
「見せる」と「測る」が同じ分解を共有し、band_snr が両方に代償を貼り付ける構図です。データ種は video → table → video / pairs / image2d で繋がります。
flowchart LR
A[カメラのクリップ または synthesize_translation] -->|video| B[temporal_band_power どこが何Hzで揺れているか]
A -->|video| C[band_snr 静止・帯域内・帯域外パワーの分解]
A -->|video| D[motion_magnify 帯域内位相を alpha 倍]
A -->|video| E[phase_displacement 画素ごとの dx dy と rank]
D -->|table| F[増幅動画 と image_snr_change_db]
D -->|table| G[motion_snr_out_db 利得補正済み・決して上がらない]
E -->|table| H[displacement_series 全画面の変位波形]
H -->|pairs| I[dsp.spectrum 共振周波数を読む]
C --> F
C --> G
増幅の内側は分解 → 位相編集 → 再合成の一本道で、alpha = 1 が厳密な恒等であることが全体の検算になります(実測 7.77e-16)。
flowchart LR
P[1 フレーム image2d] --> Q[complex_steerable_decompose 向きつき解析信号]
Q -->|table| R[帯域ごとの局所位相]
R --> S[temporal_bandpass 時間帯域で選ぶ]
S --> T[位相に alpha-1 を掛けて足し戻す]
T --> U[complex_steerable_reconstruct tight frame で厳密復元]
Q -.alpha=1 なら往復は機械精度.- U
import motionmag as M
# 見えない振動: 0.2 px を 4 Hz で、64 frame / 32 fps のクリップに入れる
clip = M.synthesize_translation((64, 64), 64, amplitude_px=0.2,
frequency_hz=4.0, fps=32.0, noise_sigma=0.01)
# (a) 測る — 何画素動いたか
series = M.displacement_series(clip, 3.0, 5.0, 32.0) # (T, 2) の dx, dy
print(abs(series[:, 0]).max()) # 0.2003278... px 真値 0.2
# (b) 見せる — 8 倍に増幅し、その代償を同じ返り値で受け取る
r = M.motion_magnify(clip, alpha=8.0, f_lo=3.0, f_hi=5.0, fps=32.0)
print(r["snr_in"]["image_snr_db"], r["snr_out"]["image_snr_db"]) # 29.25 -> 13.74 dB
print(r["motion_snr_change_db"]) # -2.18 dB 上がることは無い
print(r["band_power_ratio"]) # 0.629 線形からの外れ具合
print(abs(M.displacement_series(r["video"], 3.0, 5.0, 32.0)[:, 0]).max())
# 1.6029... px = 8 * 0.2
sum_k cos^(2K-2)(theta - pi k/K) = K*C(2K-2, K-1)/2^(2K-2)(theta に依らず一定)はここから。| 主張 | 実測 |
|---|---|
| 分解 → 再合成の往復 | 6.66e-16(64x64 既定)/ 7.22e-16(31x37、奇数・非正方)/ 最悪 7.77e-16(段数 1-8 x 向き 1-16 の全 128 通り) |
| 増幅後の変位 = alpha * d | 最悪絶対誤差 1.58e-14、最悪相対誤差 3.26e-13(alpha は 0, ±0.5, ±1, ±2, ±4, 8, 20 / d は 0.01 から 0.5 px) |
| alpha = 1 は雑音クリップでも厳密な恒等 | 7.77e-16 |
| 符号が反転していない | d(+3) + d(-3) の最大絶対値 5.34e-15 |
| 帯域外の運動は増幅されない | 増幅率 1.000000000000(alpha = 1, 4, 16, 100。alpha = 100 でも波形のずれは 8.53e-14) |
| 時間帯域通過が単一成分を復元 | 4.36e-15(DC 0.5 + 4 Hz 1.0 + 12 Hz 0.3 の混合から 4 Hz だけ) |
| 帯域パワー = 振幅²/2(Parseval) | 相対誤差 3.08e-16 |
| 開口問題 | 全画素 rank 1、dx = 0.3000000000000001、dy = 0.0(厳密) |
| テスト | tests/test_motionmag.py 124 passed / 1 skipped(skip は構成上わざと線形領域の外に置いた 1 件) |
変位計測の精度と、それが止まる場所(64x64x64 / 32 fps / 8 px 格子 / 4 Hz / 無雑音):
| 真値 d [px] | k*d [rad] | 実測 [px] | 相対誤差 | reference_coherence |
|---|---|---|---|---|
| 0.001 | 0.0008 | 0.00100000 | 8.7e-15 | 1.00000 |
| 0.010 | 0.0079 | 0.01000000 | 3.1e-15 | 0.99999 |
| 0.100 | 0.0785 | 0.10000000 | 1.8e-15 | 0.99923 |
| 0.500 | 0.3927 | 0.50000000 | 6.7e-16 | 0.98091 |
| 1.000 | 0.7854 | 1.00000000 | 4.4e-16 | 0.92582 |
| 2.000 | 1.5708 | 2.00000000 | 6.7e-16 | 0.73600 |
| 3.000 | 2.3562 | 3.00000000 | 5.9e-16 | 0.51283 |
| 3.050 | 2.3954 | 3.05000000 | 2.9e-16 | 0.50252 |
| 3.100 | 2.4347 | 1.72842712 | 4.4e-01 | 0.50761 |
| 4.000 | 3.1416 | 2.10461396 | 4.7e-01 | 0.65142 |
| 6.000 | 4.7124 | 2.35441757 | 6.1e-01 | 0.62047 |
雑音下では精度は方法ではなく雑音で決まります: d = 0.5 px で sigma = 0.001 / 0.01 / 0.05 のとき相対誤差 2.2e-04 / 1.8e-03 / 1.9e-03。
上の表は d = 3.05 px まで丸め誤差の精度で、d = 3.10 px で突然壊れます。これは調整の失敗ではなく、閉形式で位置が言えます。
時間方向の位相基準に使っているのは、その帯域の時間平均 z_mean です。振幅 A、角周波数 w の正弦的変位に対して、帯域の複素応答は z(t) = c * exp(i*phi0 - i*k*A*sin(w t)) なので、時間平均は Bessel 関数の積分表示そのものになり
z_mean = c * exp(i*phi0) * (1/T) ∫ exp(-i*k*A*sin(w t)) dt = c * exp(i*phi0) * J0(k*A)
です。J0 は k*A = 2.4048255577(第 1 零点)で符号を変えます。そこを越えた瞬間に基準の振幅がゼロを通過し、位相が pi 跳ぶので、その基準に対して測った偏差はすべて pi ずれます。波長 8 px の格子では k = 2*pi/8 なので
A_限界 = 2.4048255577 / (2*pi/8) = 3.061919 px
となり、通った 3.05 と壊れた 3.10 のちょうど間に落ちます。返り値の reference_coherence は |z_mean| / mean_t|z|、すなわちこの |J0(k*A)| を帯域で混ぜたもので、実測 1.00000(0.001 px)から 0.50252(3.05 px、限界の直前)まで単調に落ちます — 基準が退化しつつあることの実行時警告です。
その先にはもっと硬い上限があります。変位は -k*d から復元するので、|k*d| が pi に達すると答えが折り返します。波長 L の帯域では |d| < L/2(上の例で 4.0 px)で、返り値の wrap_limit_px が実測した局所周波数からこの値を返します。単一帯域の位相からこれを越える手段は、どの位相ベース手法にもありません。
alpha * d が厳密なのは、1 つのサブ帯域に動く成分が 1 つだけ入っているときです。広帯域テクスチャは 1 つの帯域に複数の空間周波数を入れ、和の位相は変位の線形関数ではありません。実測(等方性雑音を Gaussian で平滑したテクスチャ、0.2 px の運動、alpha = 3):
| 平滑 sigma [px] | 計測の相対誤差 | 増幅後の不足 |
|---|---|---|
| 1.0 | 1.6e-03 | 4.8 % |
| 1.5 | 5.7e-04 | 5.5 % |
| 3.0 | 1.0e-03 | 9.1 % |
帯域を共有する空間周波数が多いほど不足が増えます。 計測側は局所周波数を画素ごとに推定するので影響が小さく、増幅側が受け止める形です。これは位相ベース処理に内在するもので、調整で消える種類の誤差ではありません。synthesize_translation の既定波長がオクターブを分けてあるのは、この条件を合成の側で目に見えるようにしておくためです。
alpha * d)を、定量化していないバイアスと引き換えにしないためです。flow と motion の仕事です。© 2026 Kazufumi Furuse — Fullseye operator documentation. Licensed under Apache-2.0.