魚眼カメラのキャリブレーション 〜 アルゴリズム
ということで昨日の下記エントリに引き続き，アルゴリズムの詳細に入っていきます．
参照元になっている論文では，多項式近似の妥当性や半径方向・円周方向の誤差等も加味した場合の完全モデルに関しても言及されていたのですが，このあたりは複雑なのとそれほど本質的ではなさそうだったので飛ばします．
１．レンズ，正規化画像座標系と画像座標系
上記エントリの最後にも載せたんですが，論文では下記の図をベースに理論展開が進められていたので，もう一度のせておきます．
２．アルゴリズム概要
この論文の問題設定では，近似方程式にシータに関する９次方程式を採用した場合，求めるべきパラメータは９個あります．p9 = [k1, k2, mu, mv, u0, v0, k3, k4, k5]
k1 : θの１次係数
k2 : θの２次係数
mu, mv : 単位長さ当たりピクセル数
u0, v0 : 主点
k3 : θの５次係数
k4 : θの７次係数
k5 : θの９次係数
（ただし，高次係数をキャリブレーションの初期から求めようとすると自由度が高すぎてうまく行かない可能性が高いので，論文では下記ステップ①〜③において，［k1, k2, mu, mv, u0, v0］を計算し，この値を初期値として非線形最適化を実施するとき（ステップ④）において高次係数の計算を実施しています．）
アルゴリズムは大きく分けて４つのステップから構成されてます．ざっくりとは．．．．
①．内部パラメータ（カメラパラメータ・レンズパラメータ）の初期化
キャリブを開始するにあたって，ひとまずの初期値決めをします．初期値は正解に近いほうがもちろんいいですが，ここではメーカー提供の焦点距離とか画像サイズとか，ざっくりわかる値を用いて初期値をエイヤで決めます．
②．キャリブボード検知点の再投影とホモグラフィーの計算
キャリブボード上の点が写っている座標と初期値からホモグラフィー行列を計算します．ここでは，一枚のキャリブボード写真に写った７０点の投影誤差の最適化を実施する過程でパラメータ初期値の改善を行います．
③．外部パラメータ（キャリブボードとカメラの位置関係）の計算
カメラとキャリブボードの位置関係を求めます．この過程で，カメラとキャリブボードの相対位置関係 R, T (i = 1, 2, .... N) が求まります．
３．キャリブボードと工作
ちなみに，キャリブボードはみんな大好きOpenCVのチェッカボードを使う前提です．円のキャリブボードも試したんですが，どうもなかなかうまく行かず，今回は諦めました．僕が使ったキャリブボードは7x10のキャリブボードで，拡大してA3に印刷したものを，ホームセンターで売っている厚さ3mmのアクリル板に貼り付けて実施しました．
ちなみに，キャリブボードのアクリル板は現地調達です（笑）「こいつ，沖縄まで行ってなにやってんの？」と思われているでしょうが，，，そんなこと気にも止めません(ﾉ´∀｀*)
みんなの味方，サンキュー真嘉比店です（笑）．那覇に来て急遽キャリブレーションをしないといけなくなったら，ここへどうぞ！
goo.gl
それでは，，，アルゴ詳細です．．．．
４．”①．内部パラメータ（カメラパラメータ・レンズパラメータ）の初期化”
内部パラメータの初期値設定ですが，論文では下記のように設定するように書かれていました．
k1, k2 : メーカ提供のレンズモデルに対して，多項式モデル（r(θ) = k1*θ + k2*θ3）が一番フィットする k1, k2 を手計算．
mu : 単位長さあたりのピクセル数
mv : 単位長さあたりのピクセル数
u0 : 画像横サイズ / 2.0
v0 : 画像縦サイズ / 2.0
が，自分はなにも設定しなくても大丈夫でした（笑）ちなみに，OpenCVの場合，明示的に初期値を設定しなかった場合，下記のようになります．
k1, k2 : 両者とも０設定．
mu, mv : max(画像縦サイズ，画像横サイズ) / PI
u0 : 画像横サイズ / 2.0 - 0.5
v0 : 画像縦サイズ / 2.0 - 0.5
単位長さ当たりピクセル数の初期値設定がこんなに適当でもいいのかな？と思ったのですが，自分の場合は収束しました．
５．”②．キャリブボード検知点の再投影とホモグラフィーの計算”
キャリブボードの検知結果からは，一枚当たり７０点の画像座標（ximg）がわかります．この画像座標を，初期値（k1, k2, u0, v0, mu, vu）を用いてレンズ歪み無しの正規化画像座標(sXworld)に変換します．下図のようなイメージです．歪みモデルを r(θ) = k1 * θ + k2 * θ 3 と仮定してますが，こいつは歪み無し -> 歪み有りの方向に作用する関数なので，θを求めるためにθの３次関数を解く必要があります．
歪み無し正規化画像座標が求まった時点で，この画像と実際のキャリブボードが射影変換で結ばれます．実際のキャリブボードは対応点が７０点存在するので，ホモグラフィー行列を求めるために線形方程式を解いてあげます．イメージとしては下記のような感じですかね．
で，めでたしめでたし．というわけではなく，結局初期値は適当に決めただけなので，歪みをとった後の正規化画像座標と実際のキャリブボードは厳密には射影変換で結ばれていません．ただ，キャリブボードの点が７０点有り，これはホモグラフィー行列を求めるために必要な点数よりもはるかに多いため，この線形方程式は Over Determined System になります．ということで，ここで一発目の非線形最適化を実施します．ここで実施する最適化は一枚の画像に対してのみ実施します．
この最適化の誤差式が自分には若干トリッキーだったんですが，下記の２つを比較します．
１．キャリブボードの点に対して，ホモグラフィー行列を計算し投影点を計算．
２．歪み取り後のボード検知画像上の対応点．
１と２のなす角度をαiとし，この角度が小さくなるように最適化を使って（k1, k2, u0, v0, mu, vu）を調整していく．
ということで，簡単に図示してみました．
６．”③．外部パラメータ（キャリブボードとカメラの位置関係）の計算”
５までのステップで，それぞれの画像に対してホモグラフィー行列が求まりました．ここではその結果を元に，キャリブボードとカメラの位置関係を表す外部パラメータを計算します．論文にあった式の導出をそのまま書いてますが，結局キャリブボード，正規化画像座標系ともにZ方向が０となるように座標を設定しているので，Z方向の計算が消えて，ホモグラフィー行列を当てはめれば求めたい式が成立するようになります．で，Z軸に関しては，X軸とY軸が求まればそれに垂直な方向に一意に決まってしまい，なんだか騙されたような感じですが，キャリブボードと正規化画像の位置関係が求まってしまいました．
ただし，上記の計算だと誤差があるせいで求まる回転行列が完全な直交行列にならないので，さらに特異値分解を使ってできるだけ直交行列に近くなるように再計算します．（このあたりはピントきてないので，次のエントリでOpenCVの実装とともに見ます．）
７．”④．非線形最適化を用いた再投影誤差の最小化”
ふー．やっと最後まで来ました．ここまでのステップでは，ここの画像に対してのパラメータ計算・最適化でしたが，最後の総仕上げとしてすべての画像セットに対して最適なパラメータになるよう，非線形最適化をします．今までのステップで，各画像における外部パラメータ（キャリブボードとカメラの位置関係（R, t））が求まり，また内部パラメータ（[k1, k2, mu, mv, u0, v0]）が求まりました．この時点でキャリブボードの各点の三次元位置がわかるので，これを仮想的に画像に投影します．この投影座標と実際の写真に写っている座標の差の二乗和が最小になるように計算します．
ということで，明日からはOpenCVの実装とともに実装編に突入です！
（きちんとまとめないと頭に残らないし，きちんとまとめるとめっちゃ時間かかるし．．．．なかなかうまいやり方がないもんですかね．．．．）
