fullseye

実写に合わせて物理モデルの係数を決める ―― 多変量解析から入る sim-to-real 較正

English version: A Good Fit Is Not a Right Answer

利得と照明の谷の底は 1 本の線になる

↑ 照明の強さを決め打ちにして、カメラの利得だけを合わせ直した点(橙)。どの点も真値の線(破線、利得 x 照明 = 900)の上に乗ります。画像はどの点でも同じように説明できるので、画像からは「どれが正しいか」を選べません。

この記事で扱うこと

シミュレーションの絵と実物の写真が合わないとき、模型の係数(焦点距離、屈折率、膜厚、カメラの利得など)を写真に合わせて直したくなります。このシリーズの出発点は、著者が書いた次の 2 行です。

実物の画像を与えられたら、誤差をフィードバックして計算に使う係数を自動的調整するようなことができるといいな

多変量解析みたいなことから入るんじゃなかろうか?

前半はよくある「最小二乗で合わせる」です。後半が要点で、合わせる前に「どの係数が写真から決まるのか」を調べておかないと、合わせた結果を信用できません。この記事では、その理由を 3 つの具体例で示します。

先に結論を書きます。

例はすべてオープンソースの画像計測ライブラリ Fullseye(https://github.com/furuse-kazufumi/fullseye)にあるツアーで再現可能です。記事の数字は、次の 2 本を実行した印字と照合してあります。

py -3.11 examples/simcal_tour.py
py -3.11 examples/simcal_checks_tour.py

用語

用語 かみくだくと この記事での使い方
つまみ(knob) 模型の中で、写真に合わせて動かす係数 照明の強さ、利得、屈折率 n、膜厚 d、焦点距離 f
χ²/自由度 残差を雑音の大きさで割って二乗し、平均したもの 1 なら「残差は雑音そのもの」。60 なら雑音の 60 倍ずれている
交絡 2 つのつまみが写真を同じように動かすこと 決まるのは片方ずつではなく組み合わせ(積など)だけ
同定性 写真からつまみが 1 つの値に決まるかどうか 決まらないなら、値を返さずに拒否する
保留画像(holdout) 当てはめに使わず、取り置いた写真 当てはめた模型が「初めて見る写真」でも合うかを見る
残差 PCA 残差の写真を主成分に分けること 雑音より大きな「形」が残っていれば、模型に足りないものがある

例 1: 利得 x 照明は、積しか決まらない

かみくだき

写真を明るくするつまみが 2 つあるとします。照明の強さと、カメラの利得です。照明を 2 倍にして利得を半分にすると、写真はまったく同じになります。写真から決まるのは「照明 x 利得」の積だけで、照明と利得を別々には決められません。

当たり前に聞こえますが、最小二乗のソフトに 2 つとも渡すと、たいていは何かしらの値が返ってきます。その値は出発点しだいで変わる 1 点で、真値とは関係がありません。

図と数字

冒頭の図は、照明を 0.4〜4 の 13 点に決め打ちにし、そのたびに利得を合わせ直した結果です。真値は照明 1.5、利得 600(積 900)です。合わせ直した利得はどれも「利得 = 900 / 照明」の線に乗ります。

次の図は、そのときの当てはまりの良さ(χ²/自由度)です。どの点も 0.94 前後で、雑音どおりの 1 を少し下回ります。照明を 0.4 にしても 4 にしても、写真の説明力は同じです。

線の上ならどこでも同じだけ合う

Fullseye の当てはめ(simcal_fit)は、この場合に値を返さずに拒否し、決まる組み合わせだけを報告します。

決まらないつまみがあるときは、どれかを別の手段(照度計、カタログ値など)で固定するか、そのつまみを動かす撮り方に変えます。単波長で撮った膜も同じ形で、屈折率 n と膜厚 d は別々に決まらず、決まるのはほぼ n·d(光路長)だけです。

例 2: 干渉次数を取り違えた当てはめは、相対の門を全部通る

かみくだき

シリコンの上に透明な膜(n = 1.46、d = 520 nm)があり、400〜700 nm の 31 波長で反射を撮ります。反射率は波長に対して波打ち、山と谷の数は膜厚で決まります。χ² の地形には干渉の次数ごとに谷が並ぶので、出発点が悪ければ、行き先は隣の次数の谷です。

落ちた先でも、当てはめは「それなりに」合って見えます。問題は、よく使われる検証がどれも相対的な比較だということです。

図と数字

隣の干渉次数に落ちた当てはめ

小さな計画(出発点 2 つ)から始めた Levenberg-Marquardt は、n = 1.329、d = 356.5 nm に落ちました(破線)。χ²/自由度は 60.4 です。出発点を 16 個にした multistart は、n = 1.4591、d = 520.7 nm(χ²/自由度 0.93)に着きました(実線)。

この 2 つの結果を、取り置いた写真(保留画像)を使って simcal_validate の 6 つの門にかけたのが次の表です。

門 比べるもの 隣の次数 正しい次数
1. 初期値より良いか 当てはめ前の χ² 通る 通る
2. 保留 / 学習 ≤ 1 + τ 学習の写真と保留の写真の誤差の比 通る 通る
3. 利得と黒だけの模型より良いか いちばん単純な模型 通る 通る
4. 範囲・事前・端 値が範囲の端に張りついていないか 通る 通る
5. つまみ ≤ 画像の方向の数 過適合 通る 通る
6. 絶対: χ² の上側確率 p ≥ 0.001 雑音の模型 落ちる 通る

隣の次数の当てはめは、初期値より良く、保留画像でも学習と同じ程度に合い(保留の χ²/自由度は 53.1 で、学習の 60.4 より小さい)、単純な模型にも勝ちます。1〜5 番の門はどれも「何かと比べて良いか」を見ているので、全部通ります。

「雑音の何倍ずれているか」を見るのは 6 番目の門だけです。学習で χ²/自由度 60.4、保留で 53.1 は、雑音どおりなら起きる確率がほぼ 0(p < 0.001)なので、ここで初めて拒否されます。

この門には弱点もあります。雑音を大きく見積もりすぎると(分散を 4 倍にすると)p は 1 に張りつき、門は通ってしまいます。雑音の過大評価は χ²/自由度が 1 よりずっと小さいことでしか見えません。雑音の模型はカメラの測定(EMVA 1288 のような線形カメラの模型)で決めておく必要があります。

例 3: 族の外のデータは、残差 PCA が形で知らせる

かみくだき

模型の「族」とは、つまみをどう動かしても作れる絵の集まりです。光軸が中心からずれたレンズの周辺減光は、中心対称の cos⁴ の模型ではつまみをどう動かしても作れません。族の外のデータです。

それでも当てはめは数を返します。足りない分をどれかのつまみが吸い込んで、ゆがんだ値で合わせるからです。

図と数字

光軸を (5, −4) px ずらしたレンズで平らな板を 4 段の露光で撮り、中心の cos⁴ 模型(つまみ: 焦点距離 f、利得、黒レベル)で合わせました。

族の外: 光軸のずれたレンズ

左から、撮った平板、残差(雑音で割ったもの)、4 枚の残差の主成分分析の第 1 成分、真のずれ(ずれたレンズの絵 − 中心の模型の絵)です。第 1 成分の形は真のずれと相関 0.99 で一致します。

残差 PCA の特異値

雑音だけなら、特異値はランダム行列の理論で決まる端(Tracy-Widom 分布の 99 % 点、ここでは 42.6)の下に収まります。光軸がずれたレンズでは第 1 特異値が 444 で、端の 10 倍を超えています。光軸がずれていない写真(族の中)では、4 つとも端の下です。simcal_discrepancy はこの場合「雑音を越える構造が 1 成分ある。不確かさは楽観側」と注記します。

残差 PCA は「何が足りないか」の形を示します。この例なら、第 1 成分から「減光の中心がずれている」と読み取れるはずです。そこから模型に光軸の位置のつまみを足すかどうかを決めます。

どう使うか

  1. 合わせる前に、つまみごとの感度を測り(simcal_screen、Morris のふるい)、効かないつまみと交絡している組を見つける。
  2. 決まらない組があれば、撮り方を変える(波長を増やす、チャートを撮る)か、どれかを外部の値で固定する。
  3. 当てはめる(simcal_fit)。谷がいくつもある問題では multistart を使う。
  4. 取り置いた写真で 6 つの門にかける(simcal_validate)。絶対の χ² の門を外さない。
  5. 残差 PCA(simcal_discrepancy)で、模型に足りないものの形を見る。
import fullseye as fs
L = fs.ledger

obs = L.simcal_observe(model, images)          # 模型と写真を組む
screen = L.simcal_screen(obs)                  # 効く / 効かない / 組でしか効かない
fit = L.simcal_fit(obs)                        # 決まらないつまみがあれば拒否
checked = L.simcal_validate(fit, holdout)      # 6 つの門
report = L.simcal_discrepancy(checked)         # 残差 PCA
print(L.simcal_report(report)["text"])

model の組み立て方は examples/simcal_tour.py の先頭にあります。

向かないこと

第 2 回: 実写の画像で合わせる、収束まで解く模型を合わせる

この回で扱うこと

第 1 回の例はすべて合成の写真でした。真値を自分で決められるので、門が正しく鳴るかを確かめるには向いています。ただし、実写では雑音の大きさも、模型に足りないものも、こちらで選べません。この回では実写で較正した 3 つの例と、観測は合成でも順方向の模型が反復計算で解かれる 1 つの例を扱います。

先に結論を書きます。

数字は、次の 4 本の実行の印字と、図を作る tools/gen_simcal_article_figs.py が書く numbers.json(part2 の項)の両方に合わせてあります。実写の画像は、各 PoC が公開されている場所から取得してキャッシュします(薄膜のスペクトルだけはリポジトリに同梱です)。

py -3.11 examples/poc_simcal_thin_film.py
py -3.11 examples/poc_simcal_telecentric.py
py -3.11 examples/poc_simcal_pixel_shift.py
py -3.11 examples/poc_pocket_interreflection.py

用語(この回で新しく出てくるもの)

用語 かみくだくと この回での使い方
プロファイル尤度 1 つのつまみを動かさずに固定し、他を合わせ直したときの χ² の増え方 増え方が 3.84 未満の範囲が 95 % 区間
純粋誤差 模型を使わずに、近くの点どうしの食い違いだけから見積もった雑音 模型の足りなさが混ざらない雑音の物差し
ジャックナイフ データを 1 まとまりずつ抜いて解き直し、答えのばらつきから不確かさを見る方法 リングを 1 本ずつ抜く
位相相関 2 枚の画像のずれを、周波数領域の相関の山から読む別の方法 simcal の答えとの突き合わせ
相互反射 壁で反射した光が別の面を照らすこと 凹んだ部品の底の明るさの半分近く
収束解 反射を、増える分が雑音に埋もれるまで数えた結果 真値の係数で描いた観測と合う

例 4: 実測の薄膜は、狭い区間でも基準を外す

かみくだき

シリコンの上の酸化膜(SiO₂)の反射スペクトルを、顕微分光で 440〜950 nm にわたって測った実測データがあります。論文が示す基準の膜厚は 98.5 nm です。第 1 回の例 2 と同じ手順で、膜厚 d、屈折率の分散(Cauchy の A と B)、対物レンズの開口数 NA、利得の 5 つのつまみを当てはめました。

雑音の大きさは、膜のない裸のシリコンを別に測ったデータから決めました(反射率で σ = 0.00054)。当てはめる写真自身の残差から雑音を決めると、模型の足りなさまで雑音に数えてしまうからです。

図と数字

5 つのつまみは、このスペクトルからは互いに肩代わりできます(局所の共線性)。そのため simcal_fit は点推定を返さずに拒否し、つまみごとの区間で答えます。膜厚 d、A、利得は決まり、B と NA は上限だけが言えます(B は 1977 nm² 以下)。

膜厚のプロファイル尤度が次の図です。95 % 区間は [94.6, 99.0] nm(幅 4.4 nm)で、基準の 98.5 nm を含みます。取り置いた 328 波長で 6 つの門をすべて通りました(絶対の χ²/自由度は学習 0.29、取り置き 0.34)。状態は partial です。B と NA が上限しか言えないため、すべてのつまみの点推定は返さない、という状態です。

実測の膜厚のプロファイル

↑ 水色が 5 つのつまみの模型です。破線の 3.84 を越えない範囲が 95 % 区間で、点線の基準 98.5 nm がその中にあります。橙は「直入射で、屈折率は材料表の値」という貧しい模型で、区間は [93.02, 93.07] nm と細く、基準から約 5 nm 外れています。

χ²/自由度が 1 より小さい(0.29)のは、雑音の床を大きめに見積もっているためです。裸のシリコンの測定から決めた σ は、波長ごとの 2 階差分で測った白色の分(0.00014)より大きい値です。

貧しい模型は、同じスペクトルを見た目ではよく説明します。残差を雑音で割ると違いが見えます。

残差を雑音の床で割る

↑ 水色が豊かな模型、橙が直入射・材料表の n の模型の残差です。橙は全帯域でうねり、±1 σ の線(破線)から出ます。

貧しい模型を門にかけた結果は次のとおりです。

最後は、同じ膜を実物の顕微鏡写真(カメラの機種は不明、JPEG)で較正した結果です。写真から決まるのは光源の色温度と周辺減光だけで、膜厚、屈折率、利得、ガンマは組み合わせでしか決まりません。屈折率を材料表の値に固定すると膜厚は d = 86.65 nm(区間 [86.43, 86.86])と狭く出ますが、基準から 12 % 低く、χ²/自由度は 11.7、残差 PCA は 2 成分の構造を見つけて拒否します。模型には CIE 1931 の等色関数しかなく、実際のカメラの分光感度が入っていないことが原因です。

向かないこと

例 5: テレセントリックの実画像は、入射瞳を「遠い」までしか言えない

かみくだき

物体側テレセントリックのレンズは、物体の距離が変わっても像の倍率が変わりません。完全なテレセントリックでは、入射瞳が無限遠にあります。実際のレンズは入射瞳までの距離 L が有限で、その逆数 1/L(テレセン度)が 0 から離れるほど、普通のレンズに近づきます。

29.5° 傾けた板(同心円の間隔 1 mm)をテレセントリックのレンズで撮った実画像に、板の傾き、向き、中心の位置、倍率(px/mm)、リング半径の基準、半径方向の歪曲 k1、1/L の 8 つの係数を持つ撮像の模型を当てはめました。観測は「穴の中心から角度 a の測定線の上で、リング k の中心線までの距離」です。8 本のリングの点のうち、測定線 0°, 4°, … の 330 点を学習に、2°, 6°, … の点を取り置きに使いました。

図と数字

まず、雑音を模型と独立に決めました。測定線 4° おきの点に局所の 2 次式を当てはめ、その残差(純粋誤差)から σ = 1.12 px です。この雑音で合わせると、残差は 1.35 px あり、絶対の門が落ちます(χ²/自由度は学習 1.51 で p = 7.5e-9、取り置き 1.81)。simcal_validate は数を取り下げます。

リングごと・30° ごとの残差の平均 47 個を純粋誤差の標準誤差で割ったのが次の図です。2 乗和は 96 で、雑音だけなら χ²(47) に従うはずの値の p は 3.2e-5 です。リングごとの滑らかな食い違いが残っています。板の印刷や縁の読みの癖に似た形です。

リングごとの残差の平均

↑ 横軸はリング 1〜8 の順に並べた 30° ごとの区間、縦軸は残差の平均を純粋誤差の標準誤差で割った値です。破線は ±2 です。

この残差を残差 PCA にかけると white(構造なし)と出ます。リング x 測定線の行列の形では、食い違いは主成分として立ちません。ここでは絶対の門と、区間ごとの平均が止めています。

食い違いを白い分散として加え(σ = 1.37 px)当てはめ直すと、6 つの門を通ります。ただし数を出す前に、不確かさが正直かを別の方法で確かめました。リングを 1 本ずつ抜くジャックナイフです。1/L でのジャックナイフと Laplace 近似の比は 0.77 で、同じ大きさでした。

入射瞳までの距離の逆数

↑ 水色が当てはめた 1/L(線は ±2σ)、橙が焦点距離 100 mm(入射瞳 302 mm)と 25 mm の通常のレンズの値です。縦軸は対数です。

通常のレンズ(焦点距離 100 mm で 1/L = 3.31 /m、25 mm で 12.9 /m)は、この板で否定できます。一方、1/L は 0 から 2.7σ 離れていて、完全なテレセントリックとも言い切れません。光軸の位置を画像の中心と仮定していることや、リングの癖と混ざっている可能性があり、3σ に届かない差なので、どちらとも決めていません。

なお合成で同じ手順を確かめたところ、リングごとに滑らかな食い違いを入れた板では、6 つの門を通っても、ジャックナイフの不確かさが Laplace の 2.4 倍ありました(食い違いのない板では 1.0 倍)。Laplace は楽観的です。食い違いが残るときは、門を通ったことだけでは不確かさを信じられません。

向かないこと

例 6: ピクセルシフトの実データは、決まる組だけ突き合わせられる

かみくだき

ピクセルシフトは、撮像素子を動かして同じ場面を 4 枚撮り、画素の隙間を埋めて解像を上げる撮影です。動かした量(ずれ)は公称値ですが、実機では公称から外れます。そのずれを、撮れた 4 枚の画像から当てはめます。

センサのカラーフィルタは配列(RGGB)のため、枚によって同じ画素が見る色が違います。そこで、斜めの公称のずれ(±1, ±1 画素)の組を選びます。この組では緑(G)どうしが重なるので、色の違いが入りません。G だけを 45° 回した格子に並べ、1 枚目を参照にして 3 次スプラインで補間し、もう一方の枚の G の格子を予測します。格子を 16 個の升に分け、升ごとにずれ(dr, dc)を当てはめました。組は 1 枚目から 3 枚目と、2 枚目から 4 枚目の 2 つです。

図と数字

雑音は、紙の明るい所で同じ点を見る G の組の差の分散から決めました(K = 0.363 DN/電子)。

当てはめの結果は、32 個の升すべてで dr と dc が決まりました。別の方法である位相相関(副画素まで読む)と、升ごとに突き合わせます。

升ごとの差

↑ 升ごとの差(simcal の答え − 位相相関の答え)です。左の 16 点が dr、右の 16 点が dc で、橙が 1 枚目から 3 枚目、水色が 2 枚目から 4 枚目の組です。破線は差 0 です。

升ごとの差の平均は、1 枚目から 3 枚目の組で (+0.0016, −0.0011) 画素、2 枚目から 4 枚目の組で (+0.0019, −0.0007) 画素です。対応のある t 検定の p は、各組の dr と dc のうち小さい方で 0.22 と 0.30 です。0 と区別できません。個々の升の差は最大 0.036 画素まで散ります。χ²/自由度の中央値は 1.55 と 1.02 で、1 に戻りません。雑音の K を 1 点から決めており、画素ごとの感度のむらを分けていないためです。

もう 1 つの組、1 枚目から 2 枚目は、G が R と B の画素を見る組です。同じ模型(モザイク全体を参照にする)で当てはめると、χ²/自由度は 31.1 で、dr は範囲の端に張りつき、simcal は値を取り下げます(状態は partial)。同じ模型の斜めの組は 1.84 です。合成の画像で同じ組を試すと、同じ模型で値が返ります。だから違いは、模型にない物理の分です。

つまり、この無彩色のチャートからは、2 つの組の間のずれは決まりません。ずれが決まらないという結果を、数を返さないことで示せています。

向かないこと

例 7: 反射を数えない模型は、あり得ない値で合う

この例だけは、観測が合成です(実写の画像は使っていません)。それでも第 2 回に入れたのは、順方向の模型が反復計算で解かれる場合に、第 1 回と同じ門がどう鳴るかを見るためです。

かみくだき

凹んだ部品(加工ポケット、内角)では、向かいの壁で反射した光が底を照らします(相互反射)。直接光だけを数える模型には、この分がありません。足りない明るさは、反射率や照明の強さに吸収されて、もっともらしく合ってしまいます。

28 x 22 x 14 mm のブロックに 20 x 14 mm、深さ 10 mm のポケットを開けた部品(反射率 0.7)を、ドーム照明と線形カメラで撮った合成の観測に、3 つの模型を合わせました。ポケットの底の明るさの平均 48 %(ポケット内の画素で最大 62 %)が相互反射です。

つまみは反射率と照明の強さの 2 つで、真値は 0.70 と 1.00 です。

図と数字

3 つの模型の結果

直接光だけの模型は、反射率 1.32 で合います。1 を超える反射率は、入る光より出る光が多いことを意味するので、物理としてあり得ません。χ²/自由度は 23.95、残差 PCA は 4 成分の構造を見つけて、絶対の門と「利得と黒だけの模型より良いか」の門が落とします。1 回反射の模型は反射率 0.873 で、χ²/自由度 7.89 です。収束解は反射率 0.6996、照明 0.9999 と真値に戻り、χ²/自由度は学習 1.01、取り置き 1.04、残差は白色で、6 つの門をすべて通ります。

反射を何回数えれば足りるかも、雑音から決められます。

反射の回数と収束解との差

↑ 橙が、k 回で止めた像と収束解との実際の差(像の最大値に対する比、対数)、水色が反復の増分から出した GCI 型の見積もりです。破線は観測の雑音(最も明るい画素で 1.32 %)です。

雑音を許容値にすると、必要な反射の回数は 7 回です(そのときの実際の差は 0.42 %)。見積もりを信じてよいのは、隣り合う増分の比が 2 回続けて 5 % 以内でしか動かない回(漸近域)だけです。その前の k = 3 では、見積もり 0.0501 が実際の差 0.0523 を下回りました。楽観した見積もりなので、本文にも残します。

実写との接点は、数値の引用だけです。飛行時間(ToF)カメラで凹んだ 3 つの形状を 2 つの材料で測った論文の表 2 では、1 回反射を足すと平均絶対誤差が直接光だけの 0.20〜0.47 倍になります。この例の強度画像では同じ比が 0.53 でした(収束解の誤差は 22.4 DN で、雑音だけの目安 21.9 DN と並びます)。向こうは距離、こちらは強度なので、比べているのは「直接光だけでは足りない量」の向きと大きさだけです。論文の値は改変せずに引用し、比は自分で計算しました。

向かないこと

第 2 回のまとめ方

  1. 雑音は、模型とは別の測定から決める(膜のない基板、同じ点を見る組、近くの点どうしの食い違い)。
  2. 当てはめる前に、決まる量と決まらない量を分ける。決まらない量は、区間か上限で答える。
  3. 絶対の χ² の門と残差 PCA を外さない。区間が狭いことは通過の理由にならない。
  4. 食い違いが残るときは、不確かさをジャックナイフなど別の方法で確かめてから数を出す。
  5. 反復で解く模型は、回数を雑音の大きさから決め、直接光だけの模型と並べて見る。

出典(実写の画像と引用した数値)

CC BY 4.0: https://creativecommons.org/licenses/by/4.0/

シリーズの現在地

回 題材 状態
1 合成の写真で、決まらない係数・取り違えた次数・族の外を見つける 公開
2 実写の薄膜・テレセントリック・ピクセルシフトと、収束まで解く相互反射の模型(上の第 2 回の節) 追記
3 第 2 回で決めきれなかった点(色をまたぐずれ、入射瞳のテレセン度、薄い膜の分散) 予定

次回予告

第 2 回は、決まらないものを決まらないと答える所までで終わりました。残った点は 3 つです。ピクセルシフトで色をまたぐ組のずれ、テレセントリックの 1/L が 0 から 2.7σ 離れていること、薄い膜で分散を無視した模型を拒否できないことです。次回は、このうち実データで決められるものを調べます。