JP4970118B2 - カメラ校正方法、そのプログラム、記録媒体、装置 - Google Patents

Description

本発明は、カメラで撮った複数枚の画像から、カメラの位置姿勢、カメラの内部パラメタ等を推定するカメラ校正技術に関する。

カメラの光学系に歪みがない理想的な系を考えた場合、カメラ撮像素子上の点ｍ＝［ｕ，ｖ］^Ｔと空間中の点Ｍ＝［Ｘ，Ｙ，Ｚ］^Ｔとの間には、レンズの焦点距離、カメラ撮像素子の画像ピッチ(水平走査方向、垂直走査方向)、光学系の光軸と撮像素子との交点座標といったカメラ内部パラメタ、カメラの位置や姿勢といったカメラ外部パラメタを用いて、

ここで、α_Ｃは、カメラ撮像素子の水平走査方向の画素ピッチと、カメラレンズの焦点距離の、積の逆数、β_Ｃは、カメラ撮像素子の垂直走査方向の画素ピッチと、レンズの焦点距離の、積の逆数を、水平走査方向と垂直走査方向の成す角の正弦（ｓｉｎ）で割った値、γ_Ｃは、カメラ撮像素子の垂直方向の画素ピッチと、レンズ焦点距離の、積の逆数を、水平走査方向と垂直走査方向の成す角の正接（ｔａｎ）で割った値、ｕ_Ｃ，０は、光学系の光軸と撮像素子との交点のカメラ撮像素子の水平走査方向の座標、ｖ_Ｃ，０は、光学系の光軸と撮像素子との交点のカメラ撮像素子の垂直走査方向の座標、Ｒ_Ｃはカメラ回転行列、ｔ_Ｃはカメラ平行並進ベクトル、のような関係が成り立つ。なお、上記の文章中に記載したｍ＝［ｕ，ｖ］^Ｔ、Ｍ＝［Ｘ，Ｙ，Ｚ］^Ｔの左辺の「ｍ」、「Ｍ」は本来は太字で記載すべきものであるが、文章中では太字が使えないため通常の文字で「ｍ」、「Ｍ」と記載した。以下の文章中の式で用いる文字についても同様である。

ところが実際には、光学系には一般に「樽型歪み」や「糸巻き型歪み」と呼ばれる歪曲歪みが存在しており、空間中の点Ｍの理想的なカメラ撮像素子上の射影点ｍと、実際にカメラ撮像素子上に結像される座標（すなわち、レンズで歪んだ場合に結像される座標）

とは一致しない。なお、数２のイメージデータに示した記号を文章中では「ブレーヴェｍ」と記載する。

また、一般にカメラ内部パラメタ、カメラ外部パラメタも厳密な値は既知ではないので、理想的なカメラ撮像素子上の射影点ｍから実際にカメラ撮像素子上に結像される座標ブレーヴェｍへの変換関数Ｄ［．］とともに、カメラ内部パラメタ、カメラ外部パラメタも、空間中の観測点の座標と撮像素子上の結像点の座標の組から推定する必要がある。

このような歪曲歪みのパラメタも含めたカメラパラメタの推定方法として、従来、以下のような方法が行われていた（たとえば、非特許文献１）。

図１は空間中にある平面上の観測点とカメラ撮像素子上の結像点の関係を示す図であり、１０１は空間中の平面Ｂ１、１０２は平面上の観測点ｐ_Ｂ１，１、１１１はカメラＣ、１１２はレンズ、１１３はレンズ中心、１１４は撮像素子、１１５は撮像素子原点、１１６は観測点ｐ_Ｂ１，１の結像点ｐ_{Ｃ，Ｂ１，１}、を示す。

カメラパラメタの推定は、たとえば、下記のステップで行われる。
Ｓｔｅｐ１．レンズ歪みが無いとして、カメラ内部／外部パラメタを概算。
Ｓｔｅｐ１−１．空間中の平面Ｂ１とカメラ撮像素子のホモグラフィ行列を概算。
Ｓｔｅｐ１−２．カメラ内部パラメタを概算。
Ｓｔｅｐ１−３．カメラ外部パラメタを概算。
Ｓｔｅｐ２．レンズ歪み係数を概算。
Ｓｔｅｐ３．カメラ内部／外部パラメタ、レンズ歪み係数を最適化手法により高精度化。

以下、詳細に説明する。

［Ｓｔｅｐ１．レンズ歪みが無いとして、カメラ内部／外部パラメタを概算］
［Ｓｔｅｐ１−１．空間中の平面Ｂ１とカメラ撮像素子のホモグラフィ行列を概算］
平面Ｂ１上の観測点ｐ_Ｂ１，１の座標を世界座標系で表し、Ｘ_{ｐＢ１，１}＝［Ｘ_{ｐＢ１，１}，Ｙ_{ｐＢ１，１}，Ｚ_{ｐＢ１，１}］^Ｔのようにデカルト座標系で表記する。また、撮像素子上にある観測点ｐ_Ｂ１，１のカメラＣの撮像素子上の結像点ｐ_{Ｃ，Ｂ１，１}の座標を撮像素子のディジタル画像座標系で表し、ｍ_{ｐＣ，Ｂ１，１}＝［ｕ_{ｐＣ，Ｂ１，１}，ｖ_{ｐＣ，Ｂ１，１}］^Ｔのように表記する。

１０１の平面Ｂ１を世界座標系の座標［Ｘ，Ｙ，Ｚ］^ＴのＺ＝０の平面とし、その座標を２次元座標Ｘ_{ｐＢ１，１}’＝［Ｘ_{ｐＢ１，１}，Ｙ_{ｐＢ１，１}］^Ｔとしたとき、平面Ｂ１上の観測点ｐ_Ｂ１，１と結像点ｐ_{Ｃ，Ｂ１，１}のそれぞれの座標は、平面Ｂ１と撮像素子の平面で関連づけられる１つのホモグラフィ行列Ｈ_Ｃ，Ｂ１＝［ｈ_{Ｃ，Ｂ１，１} ｈ_{Ｃ，Ｂ１，２} ｈ_{Ｃ，Ｂ１，３}］を使って次のように関係づけられる。ただし、ｈ_{Ｃ，Ｂ１，１}、ｈ_{Ｃ，Ｂ１，２}、ｈ_{Ｃ，Ｂ１，３}はそれぞれ３行１列の列ベクトルである。

ここで、ベクトルの上に付いている「〜」は同次座標であることを示し、チルダｍ_{ｐＣ，Ｂ１，１}＝［ｕ_{ｐＣ，Ｂ１，１}，ｖ_{ｐＣ，Ｂ１，１}，１］^Ｔ、チルダＸ_{ｐＢ１，１}’＝［Ｘ_{ｐＢ１，１}，Ｙ_{ｐＢ１，１}，１］^Ｔであり、ｓはスカラ値である。なお、数３のイメージに記載されているような「〜」が上についた「ｍ」等を文章中では「チルダｍ」等と記載する。

平面Ｂ１上の点の座標とそれに対応する撮像素子上の結像点の座標を（３．１）式に代入して方程式をつくり、ホモグラフィ行列Ｈ_Ｃ，Ｂ１＝［ｈ_{Ｃ，Ｂ１，１} ｈ_{Ｃ，Ｂ１，２} ｈ_{Ｃ，Ｂ１，３}］を算出する。この方法としては、たとえば、下記の方法がある。まず、ｎ組の観測点の座標と結像点の座標の組を使って、次の行列Ｌ_Ｃ，Ｂ１を作る。

ここで、０＝［０，０，０］^Ｔである。

次にＬ_Ｃ，Ｂ１ ^ＴＬ_Ｃ，Ｂ１の固有ベクトルを求め、その中の最も固有値が小さい固有ベクトルをｘ_Ｃ，Ｂ１＝［ｈ_{Ｃ，Ｂ１，１} ^Ｔ ｈ_{Ｃ，Ｂ１，２} ^Ｔ ｈ_{Ｃ，Ｂ１，３} ^Ｔ］とする。最後にこのベクトルｘ_Ｃ，Ｂ１の各エレメントを並び直すことによりホモグラフィ行列Ｈ_Ｃ，Ｂ１＝［ｈ_{Ｃ，Ｂ１，１} ｈ_{Ｃ，Ｂ１，２} ｈ_{Ｃ，Ｂ１，３}］を得る。

［Ｓｔｅｐ１−２．カメラ内部パラメタを概算］
次に、ホモグラフィ行列Ｈ_Ｃ，Ｂ１＝［ｈ_{Ｃ，Ｂ１，１} ｈ_{Ｃ，Ｂ１，２} ｈ_{Ｃ，Ｂ１，３}］のエレメントを次のように表記する。

ただし、ｉは１≦ｉ≦３の整数である。

上記表記を使って、次のようなベクトルを作る。

ただし、ｊは１≦ｊ≦３の整数である。

平面の位置をｎ_２回変えてホモグラフィ行列Ｈ_Ｃ，Ｂ１＝［ｈ_{Ｃ，Ｂ１，１} ｈ_{Ｃ，Ｂ１，２} ｈ_{Ｃ，Ｂ１，３}］を作り、それらから得たｎ_２個のベクトルＶ_{Ｃ，Ｂ１，ｉ，ｊ}〜Ｖ_{Ｃ，Ｂｎ２，ｉ，ｊ}を使って、次の行列Ｖ_Ｃを作る。

次にＶ_Ｃ ^ＴＶ_Ｃの固有ベクトルを求め、その中の最も固有値が小さい固有ベクトルをｂ_Ｃ＝［Ｂ_{Ｃ，１，１} Ｂ_{Ｃ，１，２} Ｂ_{Ｃ，２，２} Ｂ_{Ｃ，１，３} Ｂ_{Ｃ，２，３} Ｂ_{Ｃ，３，３}］とする。

このｂ_Ｃを用いてカメラ内部行列Ａ_Ｃは以下のように概算される。

ただし、

ここで、α_Ｃは、カメラ撮像素子の水平走査方向の画素ピッチと、レンズの焦点距離の、積の逆数、β_Ｃは、カメラ撮像素子の垂直走査方向の画素ピッチと、レンズの焦点距離の、積の逆数を、水平走査方向と垂直走査方向の成す角の正弦（ｓｉｎ）で割ったもの、γ_Ｃは、カメラ撮像素子の水平走査方向の画素ピッチと、レンズの焦点距離の、積の逆数を、水平走査方向と垂直走査方向の成す角の正接（ｔａｎ）で割ったもの、ｕ_Ｃ，０は、光学系の光軸と撮像素子との交点のカメラ撮像素子の水平走査方向の座標、ｖ_Ｃ，０は、光学系の光軸と撮像素子との交点のカメラ撮像素子の垂直走査方向の座標、である。

［Ｓｔｅｐ１−３．カメラ外部パラメタを概算］
一方、カメラ外部行列Ｇ_Ｃ，Ｂ１をカメラ外部パラメタをエレメントに持つカメラ回転行列Ｒ_Ｃ，Ｂ１＝［ｒ_{Ｃ，Ｂ１，１} ｒ_{Ｃ，Ｂ１，２} ｒ_{Ｃ，Ｂ１，３}］とカメラ平行並進ベクトルｔ_Ｃ，Ｂ１からＧ_Ｃ，Ｂ１＝［Ｒ_Ｃ，Ｂ１ ｔ_Ｃ，Ｂ１］と定義すると、カメラ外部行列Ｇ_Ｃ，Ｂ１はこれまでの結果を用いて、

と概算される。

［ｓｔｅｐ２．レンズ歪み係数を概算］
次に、レンズ歪みを反映するためのレンズ歪み係数を計算する。従来の方法においては、レンズで歪んだ場合に結像される座標ブレーヴェｘ_{ｐＣ，Ｂ１，１}＝［ブレーヴェｘ_{ｐＣ，Ｂ１，１}，ブレーヴェｙ_{ｐＣ，Ｂ１，１}］^Ｔは、レンズ歪みがない場合に結像される座標ｘ_{ｐＣ，Ｂ１，１}＝［ｘ_{ｐＣ，Ｂ１，１}，ｙ_{ｐＣ，Ｂ１，１}］^Ｔ、レンズ歪み係数ｋ_Ｃ，１、ｋ_Ｃ，２を用いて、

としていた。ただし、本式はカメラＣのカメラ座標系で記述している。なお、数１０にイメージデータで示された（３．１０）式の左辺のｘ、ｙを文章中では「ブレーヴェｘ」、「ブレーヴェｙ」と記載する。

この変換を、レンズ歪みがない場合に結像される座標ｘ_{ｐＣ，Ｂ１，１}＝［ｘ_{ｐＣ，Ｂ１，１}，ｙ_{ｐＣ，Ｂ１，１}］^Ｔを入力とし、レンズで歪んだ場合に結像される座標ブレーヴェｘ_{ｐＣ，Ｂ１，１}＝［ブレーヴェｘ_{ｐＣ，Ｂ１，１}，ブレーヴェｙ_{ｐＣ，Ｂ１，１}］^Ｔを出力するレンズ歪み関数Ｄ［ｋ_Ｃ，１，ｋ_Ｃ，２，ｘ_{ｐＣ，Ｂ１，１}］を導入し、ブレーヴェｘ_{ｐＣ，Ｂ１，１}＝Ｄ［ｋ_Ｃ，１，ｋ_Ｃ，２，ｘ_{ｐＣ，Ｂ１，１}］と表す。

上記の式をカメラ座標系ｘ_Ｃ＝［ｘ_Ｃ，ｙ_Ｃ］とディジタル座標系ｍ_Ｃ＝Ａ_Ｃｘ_Ｃの関係を使い、レンズ歪み係数ｋ_Ｃ，１、ｋ_Ｃ，２を以下のように算出する。ただし、Ａ_Ｃのエレメントの１つであるγ_Ｃは通常非常に小さいので、計算量と精度との兼ね合いを考慮して０として計算する。

ただし、

なお、上記の式では、一つの平面上の観測点と撮像素子上の結像点の組が、１つの平面あたりｎ_１組あり、その平面をｎ_２回場所を変えてカメラで撮像した、観測点と結像点の組を全て使用している。

［Ｓｔｅｐ３．カメラ内部／外部パラメタ、レンズ歪み係数を最適化手法により高精度化］
次に、下記の式が最も小さくなるように、カメラ内部行列Ａ_Ｃ、カメラ回転行列Ｒ_Ｃ，Ｂｊ、カメラ平行並進ベクトルｔ_Ｃ，Ｂｊ、レンズ歪み係数ｋ_Ｃ，１、ｋ_Ｃ，２を最適化し、カメラ内部／外部パラメタおよびレンズ歪み係数を高精度化する。

ここで、ｎ_１は１枚の平面上にある観測点と撮像素子上の結像点の組の数、ｎ_２は平面の場所を変えて撮像した回数、Ｘ_{ｐＢｊ，ｉ}’は観測点の座標、ブレーヴェｍ_{ｐＣ，Ｂｊ．ｉ}は結像点の実測した座標、ｍ［．］は観測点Ｘ_{ｐＢｊ，ｉ}’をＡ_Ｃ、ｋ_Ｃ，１、ｋ_Ｃ，２、Ｒ_Ｃ，Ｂｊ、ｔ_Ｃ，Ｂｊを用いてカメラの撮像素子上に結像した結像点の座標を計算する投射関数であり、これまで求めてきた式を用いて、

と表される。ここでチルダｘ_{ｐＣ，Ｂｊ，ｉ}は、

によって計算される。（３．１２）式の最適化は、たとえば、Levenberg-Marquardt法により行う。

以上の手順によって、カメラ内部／外部パラメタおよび歪み係数が求められ、これを用いて歪みの補正などの各種画像処理の際のカメラ校正が行われる。

Z. Zhang,"A Flexible New Technique for Camera Calibration,"Technical Report MSR-TR-98-71, Microsoft Research, Dec. 1998. Available together with the software at http://research.microsoft.com/~zhang/Calib/.

しかしながら、従来の方法で求められたカメラ内部／外部パラメタおよび歪み係数を用いて実際に画像処理を行うと、レンズの歪みが残っている場合があった。特に近年多く使われるようになった、携帯電話のカメラなどの安価な光学系を用いたカメラで顕著であった。

上記に示すような従来のレンズ歪み補正方法は、放射方向のレンズ歪みはレンズの偏心がない場合の歪曲歪みに対応しているだけであり、レンズの光軸からの偏心がある場合の歪曲歪みについては考慮されておらず、偏心に依存した歪み量算出のためのカメラパラメタとして何を導入すればよいかの指針も明らかになっていなかったことに起因する。このため従来法では、一般のレンズが偏心しているカメラを使用した場合には歪み補正が十分なされない、という欠点があった。

この点について、図２を用いて説明する。レンズの偏心とは、光学系の基準軸に対してレンズが平行移動したり傾いたりすることにより、レンズ面の曲率中心が基準軸上からずれることを言う。図２（ａ）はこの様子を示したものであり、上図はずれ量を平行移動距離で表したものであり、下図は角度で表したものである。軸のずれがわずかな場合は、Ｅ＝ｌεであるので、平行移動Ｅと角度εは比例関係となるので、歪み係数を求める際には本質的な違いがなく、どちらで表しても良い。

図２（ｂ）は、偏心量をＥで表し、偏心した場合の歪み量について説明する図である。偏心した場合の歪み量は、無歪み時の結像位置ｌに加えて、偏心に依存しない場合の放射方向歪み（Ｖ_３‖ｌ‖^２＋Ｖ_５‖ｌ‖^４）と、偏心に依存する放射方向歪み（Ｖ_３Ｅ＋Ｖ_５Ｅ‖ｌ‖^２）（Ｅ・ｌ）ｌがある。従来は、前者の歪みのみを考慮していたが、後者の歪みもあるため、従来の方法では偏心がある場合は後者の歪みの補正が十分なされない。この理由としては、後者は（Ｅ・ｌ）の項を持つため基準軸からの方向によってその歪みの大きさが変わるという性質を持つが、従来考慮していた前者の歪みは基準軸からの方向により歪みの大きさが変わらないため、前者の歪みを考慮するだけでは、後者の歪みは十分補正できない。

基準軸からの方向により歪む量が異なることについて、結像する像の形を用いて説明する。図２（ｃ）にレンズの偏心がある場合と無い場合とでの歪み画像の違いを示す。点線は無歪み時の画像であり、一点鎖線はレンズの偏心が無い場合の歪み画像であり、実線はレンズの偏心がある場合の歪み画像である。偏心がなければ、光学系の基準軸を中心にどの角度にも同じ歪みが出て、一点鎖線のような画像となる。偏心があった場合には、矢印に示す偏心の方向に像が伸縮する。従来は、一点鎖線で示した偏心が無い場合の方向性のない歪みのみを補正対象としているが、レンズの偏心がある一般的な場合では、レンズの偏心による実線で示したレンズの偏心方向に依存した歪みが出ており、この歪みは従来方法では十分に補正することはできない。

本発明は、上記に鑑みてなされたものであり、その目的とするところは、偏心しているレンズを有するカメラにおいても、歪み補正を精度よく行うカメラ校正方法を提供することにある。

本願において開示される発明のうち、代表的なものの概要を簡単に説明すれば、以下の通りである。

第１の発明は、カメラで撮影された画像から、カメラ光学系の歪みを補正するカメラ校正装置におけるカメラ校正方法であって、校正の際に使われるレンズ歪み関数が、偏心に依存するレンズ歪みをパラメタとして有し、撮像面上のレンズの光軸中心からの距離の３乗と撮像面の各軸方向への偏心量に比例したレンズ歪みと、撮像面上のレンズの光軸中心からの距離の５乗と撮像面の各軸方向への偏心量に比例したレンズ歪みとの和を含む式として表されること、を特徴とする。

第２の発明は、カメラで撮影された画像から、カメラ光学系の歪みを補正するカメラ校正装置におけるカメラ校正方法であって、校正の際に使われるレンズ歪み関数が、偏心に依存するレンズ歪みをパラメタとして有し、一つの平面上にある観察点の座標と撮像面上の結像点の座標の組からレンズ歪みが無いとして、カメラ内部／外部パラメタを概算する第１のステップと、光学系の偏心に依存しない、撮像面上のレンズの光軸中心からの方向によらず、撮像面上のレンズの光軸中心からの距離のみにより歪み量が異なる歪みと、光学系の偏心に依存し、撮像面上のレンズの光軸中心からの方向と距離の両方により歪み量が異なる歪みの両方を足し合わせたレンズ歪み関数を用いて、３次元空間中の点のカメラ撮像面上の観測位置と、前記レンズ歪み関数を用いて計算したカメラ撮像上の位置との差の二乗が最も小さくなる最小二乗法で歪み係数を概算する第２のステップと、前記歪み関数と、カメラ位置と３次元上の観測点で規定されるエピポーラ幾何学の拘束条件を用いて、３次元空間中の点のカメラ撮像面上の観測位置と、前記レンズ歪み関数を用いて計算したカメラ撮像面上の位置との差が極小となるような最適化手法を用いてカメラ内部／外部パラメタとレンズ歪み係数を高精度化する第３のステップ、を持つことを特徴とする。

第３の発明は、第１または第２の発明のカメラ校正方法をコンピュータに実行させるためのプログラムである。

第４の発明は、第１または第２の発明のカメラ校正方法をコンピュータに実行させるためのプログラムを格納した記憶媒体である。

第５の発明は、第１または第２の発明のカメラ校正方法を実行するカメラ校正装置である。

以下、本発明の理解を助けるために、より具体的に説明するが、これにより本発明が限定されるものではない。

本発明により、レンズの偏心に依存しない従来の歪みに加えて、レンズの偏心に依存する歪みも加えたレンズ歪みの補正を行うため、高精度なレンズ歪みの補正が可能となる。

以下、図面を用いて本発明の実施の形態を説明する。

図１は空間中にある平面上の観測点とカメラ撮像素子上の結像点の関係を示す図であり、１０１は空間中の平面Ｂ１、１０２は平面上の観測点ｐ_Ｂ１，１、１１１はカメラＣ、１１２はレンズ、１１３はレンズ中心、１１４は撮像素子、１１５は撮像素子原点、１１６は観測点ｐ_Ｂ１，１の結像点ｐ_{Ｃ，Ｂ１，１}、を示す。

図３に本実施形態のカメラ校正方法の処理のフローを示す。なお、本明細書においては、カメラ内部／外部パラメタおよびレンズ歪み係数等のパラメタの算出方法をカメラ校正方法と呼ぶ。

カメラ校正に用いられるカメラ内部／外部パラメタおよびレンズ歪み係数は、以下のように求められる。以下の処理はカメラ校正装置（図示していない）によって実行することができる。カメラ校正装置は記憶装置に記憶されたプログラムとコンピュータによって構成することができ、また、そのプログラムの一部または全部に代えてハードウェアで構成することもできる。
Ｓｔｅｐ０．観測点座標と結像点座標の組の入力。
Ｓｔｅｐ１．レンズ歪みが無いとして、カメラ内部／外部パラメタを概算。
Ｓｔｅｐ１−１．空間中の平面Ｂ１とカメラ撮像素子のホモグラフィ行列を概算。
Ｓｔｅｐ１−２．カメラ内部パラメタを概算。
Ｓｔｅｐ１−３．カメラ外部パラメタを概算。
Ｓｔｅｐ２．レンズ歪み係数を概算。
Ｓｔｅｐ３．カメラ内部／外部パラメタ、レンズ歪み係数を最適化手法により高精度化。

処理の流れは従来と同じであるが、Ｓｔｅｐ２以降で求められるレンズ歪み係数の中に、レンズの偏心に依存した成分が新たに導入されている。

カメラ内部／外部パラメタおよびレンズ歪み係数（歪曲歪みの係数およびレンズ偏心の歪み係数）を求めるために、一平面上にあるｎ_１個（ｎ_１は８以上の整数）の観察点を、カメラと平面の相対位置を変えてｎ_２回（ｎ_２は３以上の整数）撮影したときの、空間中にある観察点の座標とカメラ撮像素子上の結像点の座標の組を入力として用いる。

以下、カメラ内部／外部パラメタおよびレンズ歪み係数の求め方を詳細に説明する。
［Ｓｔｅｐ０．観測点座標と結像点座標の組の入力］
平面Ｂ１上の観測点の座標と、それが撮像素子上で結像している点の座標の組を入力する。このとき、平面Ｂ１上の観測点の座標は、世界座標系で表し、また、撮像素子上で結像している点の座標は、カメラ座標軸で表す。

［Ｓｔｅｐ１．レンズ歪みが無いとして、カメラ内部／外部パラメタを概算］
［Ｓｔｅｐ１−１．空間中の平面Ｂ１とカメラ撮像素子のホモグラフィ行列を概算］
平面Ｂ１上の観測点ｐ_Ｂ１，１の座標を世界座標系で表し、Ｘ_{ｐＢ１，１}＝［Ｘ_{ｐＢ１，１}，Ｙ_{ｐＢ１，１}，Ｚ_{ｐＢ１，１}］^Ｔのようにデカルト座標系で表記する。また、撮像素子上にある観測点ｐ_Ｂ１，１のカメラＣの撮像素子上の結像点ｐ_{Ｃ，Ｂ１，１}の座標を撮像素子のディジタル画像座標系で表し、ｍ_{ｐＣ，Ｂ１，１}＝［ｕ_{ｐＣ，Ｂ１，１}，ｖ_{ｐＣ，Ｂ１，１}］^Ｔのように表記する。

１０１の平面Ｂ１を世界座標系の座標［Ｘ，Ｙ，Ｚ］^ＴのＺ＝０の平面とし、その座標を２次元座標Ｘ_{ｐＢ１，１}’＝［Ｘ_{ｐＢ１，１}，Ｙ_{ｐＢ１，１}］^Ｔとしたとき、平面Ｂ１上の観測点ｐ_Ｂ１，１と結像点ｐ_{Ｃ，Ｂ１，１}のそれぞれの座標は、平面Ｂ１と撮像素子の平面で関連づけられる１つのホモグラフィ行列Ｈ_Ｃ，Ｂ１＝［ｈ_{Ｃ，Ｂ１，１} ｈ_{Ｃ，Ｂ１，２} ｈ_{Ｃ，Ｂ１，３}］を使って次のように関係づけられる。ただし、ｈ_{Ｃ，Ｂ１，１}、ｈ_{Ｃ，Ｂ１，２}、ｈ_{Ｃ，Ｂ１，３}はそれぞれ３行１列の列ベクトルである。

ここで、ベクトルの上に付いている「〜」は同次座標であることを示し、チルダｍ_{ｐＣ，Ｂ１，１}＝［ｕ_{ｐＣ，Ｂ１，１}，ｖ_{ｐＣ，Ｂ１，１}，１］^Ｔ、チルダＸ_{ｐＢ１，１}’＝［Ｘ_{ｐＢ１，１}，Ｙ_{ｐＢ１，１}，１］^Ｔであり、ｓはスカラ値である。

平面Ｂ１の点の座標とそれに対応する撮像素子上の結像点の座標を（８．１）式に代入して方程式をつくり、ホモグラフィ行列Ｈ_Ｃ，Ｂ１＝［ｈ_{Ｃ，Ｂ１，１} ｈ_{Ｃ，Ｂ１，２} ｈ_{Ｃ，Ｂ１，３}］を算出する。この方法としては、たとえば、下記の方法がある。まず、ｎ_１組（ただしｎ_１は８以上の整数）の観測点の座標と結像点の座標の組を使って、次の行列Ｌ_Ｃ，Ｂ１を作る。

ここで、０＝［０，０，０］^Ｔである。

次にＬ_Ｃ，Ｂ１ ^ＴＬ_Ｃ，Ｂ１の固有ベクトルを求め、その中の最も固有値が小さい固有ベクトルをｘ_Ｃ，Ｂ１＝［ｈ_{Ｃ，Ｂ１，１} ^Ｔ ｈ_{Ｃ，Ｂ１，２} ^Ｔ ｈ_{Ｃ，Ｂ１，３} ^Ｔ］とする。最後にこのベクトルｘ_Ｃ，Ｂ１の各エレメントを並び直すことによりホモグラフィ行列Ｈ_Ｃ，Ｂ１＝［ｈ_{Ｃ，Ｂ１，１} ｈ_{Ｃ，Ｂ１，２} ｈ_{Ｃ，Ｂ１，３}］を得る。

［Ｓｔｅｐ１−２．カメラ内部パラメタを概算］
次に、ホモグラフィ行列Ｈ_Ｃ，Ｂ１＝［ｈ_{Ｃ，Ｂ１，１} ｈ_{Ｃ，Ｂ１，２} ｈ_{Ｃ，Ｂ１，３}］のエレメントを次のように表記する。

ただしｉは１≦ｉ≦３の整数である。

上記表記を使って、次のようなベクトルを作る。

ただしｊは１≦ｊ≦３の整数である。

平面の位置をｎ_２回（ただしｎ_２は３以上の整数）変えてホモグラフィ行列Ｈ_Ｃ，Ｂ１＝［ｈ_{Ｃ，Ｂ１，１} ｈ_{Ｃ，Ｂ１，２} ｈ_{Ｃ，Ｂ１，３}］を作り、それらから得たｎ_２個のベクトルｖ_{Ｃ，Ｂ１，ｉ，ｊ}〜ｖ_{Ｃ，Ｂｍ，ｉ，ｊ}を使って、次の行列Ｖ_Ｃを作る。

次にＶ_Ｃ ^ＴＶ_Ｃの固有ベクトルを求め、その中の最も固有値が小さい固有ベクトルをｂ_Ｃ＝［Ｂ_{Ｃ，１，１} Ｂ_{Ｃ，１，2} Ｂ_Ｃ，2，2 Ｂ_{Ｃ，１，3} Ｂ_Ｃ，2，3 Ｂ_Ｃ，3，3］とする。

このｂ_Ｃを用いてカメラ内部行列Ａ_Ｃは以下のように概算される。

ただし、

ここで、α_Ｃは、カメラ撮像素子の水平走査方向の画素ピッチと、レンズの焦点距離の、積の逆数、β_Ｃは、カメラ撮像素子の垂直走査方向の画素ピッチと、レンズの焦点距離の、積の逆数を、水平走査方向と垂直走査方向の成す角の正弦（ｓｉｎ）で割ったもの、γ_Ｃは、カメラ撮像素子の水平走査方向の画素ピッチと、レンズの焦点距離の、積の逆数を、水平走査方向と垂直走査方向の成す角の正接（ｔａｎ）で割ったもの、ｕ_Ｃ，０は、光学系の光軸と撮像素子との交点のカメラ撮像素子の水平走査方向の座標、ｖ_Ｃ，０は、光学系の光軸と撮像素子との交点のカメラ撮像素子の垂直走査方向の座標、である。

［Ｓｔｅｐ１−３．カメラ外部パラメタを概算］
一方、カメラ外部行列Ｇ_Ｃ，Ｂ１を、カメラ外部パラメタをエレメントに持つカメラ回転行列Ｒ_Ｃ，Ｂ１とカメラ平行並進ベクトルｔ_Ｃ，Ｂ１からＧ_Ｃ，Ｂ１＝［Ｒ_Ｃ，Ｂ１ｔ_Ｃ，Ｂ１］と定義すると、カメラ外部行列Ｇ_Ｃ，Ｂ１はこれまでの結果を用いて

ただし、

と概算される。

［Ｓｔｅｐ２．レンズ歪み係数を概算］
次に、レンズ歪みを反映するためのレンズ歪み係数を計算する。ここで考慮する歪みは、レンズ歪みの中で一般に歪み量の大きい、放射方向歪み（歪曲歪み）、および、レンズの曲率中心がレンズ基準軸からずれている状態、つまり、偏心している状態における歪みである。

図２にレンズの偏心と歪みの関係を示す。偏心には図２（ａ）に示すように、２種類の表示方法があり、１つは図２（ａ）の上図に示すように並進移動Ｅによる表記、もう１つは図２（ａ）の下図に示すように角度εによる表記である。ずれ量が小さい場合は、曲率半径Ｒを用いてＥ＝Ｒεとなり、平行移動Ｅと角度εは比例関係となるので、歪み係数を求める際にはどちらを用いても本質的には違いがなく、どちらを用いても良い。一般に、レンズの偏心による歪みは、平行移動Ｅと角度εを比例関係として扱ってよい程度である。

図２（ｂ）はＥを用いて歪みを示した図である。無歪み時の結像位置ｌに対して、
偏心に依存しない歪み（Ｖ_３‖ｌ‖^２＋Ｖ_５‖ｌ‖^４）と、それに加えて偏心に依存する歪み（Ｖ_３Ｅ＋Ｖ_５Ｅ‖ｌ‖^２）（Ｅ・ｌ）ｌがあり、それらを加えた

が全歪みとなる。

これらをカメラＣのカメラ座標系の各エレメントで書き下したものが下記式（８.１０）である。ただし、係数Ｖ_３、Ｖ_５、Ｖ_３Ｅ、Ｖ_５Ｅ、Ｅ＝［Ｅ_ｘ Ｅ_ｙ］はそのまま使わず、後述する（８．１１）式によって行列演算しやすいように、各項に乗じる係数を適宜組み合わせ、それらをｋ_Ｃ，１、ｋ_Ｃ，２、ｑ_Ｃ，１、ｑ_Ｃ，２、ｑ_Ｃ，３、ｑ_Ｃ，４という係数に置き換えた。

ただし、

である。

ここで、ブレーヴェｘ_{ｐＣ，Ｂ１，１}＝［ブレーヴェｘ_{ｐＣ，Ｂ１，１}，ブレーヴェｙ_{ｐＣ，Ｂ１，１}］^Ｔはレンズで歪んだ場合に結像される座標、ｘ_{ｐＣ，Ｂ１，１}＝［ｘ_{ｐＣ，Ｂ１，１}，ｙ_{ｐＣ，Ｂ１，１}］^Ｔはレンズ歪みがない場合に結像される座標、ｋ_Ｃ，１はレンズの偏心に依存しない放射方向のレンズ歪み係数で、レンズ中心からの距離の３乗に比例した歪みの係数、ｋ_Ｃ，２はレンズの偏心に依存しない放射方向のレンズ歪み係数で、レンズ中心からの距離の５乗に比例した歪みの係数、ｑ_Ｃ，１はレンズの偏心に依存した放射方向のレンズ歪み係数で、レンズ中心からの距離の３乗と、撮像素子の水平走査方向の偏心量に依存した歪みの係数、ｑ_Ｃ，２はレンズの偏心に依存した放射方向のレンズ歪み係数で、レンズ中心からの距離の５乗と、撮像素子の水平走査方向の偏心量に依存した歪みの係数、ｑ_Ｃ，３はレンズの偏心に依存した放射方向のレンズ歪み係数で、レンズ中心からの距離の３乗と、撮像素子の垂直走査方向の偏心量に依存した歪みの係数、ｑ_Ｃ，４はレンズの偏心に依存した放射方向のレンズ歪み係数で、レンズ中心からの距離の５乗と、撮像素子の垂直走査方向の偏心量に依存した歪みの係数、である。

（８−１０）式の上の２つの式のそれぞれの第２項目は偏心に関係なく現れる歪み量を表しており、第３項目、第４項目は偏心に依存する歪み量を表している。

上記の式をカメラ座標系ｘ_Ｃ＝［ｘ_Ｃ，ｙ_Ｃ］とディジタル座標系ｍ_Ｃ＝Ａ_Ｃｘ_Ｃの関係を使い、レンズ歪み係数ｋ_Ｃ，１、ｋ_Ｃ，２、ｑ_Ｃ，１、ｑ_Ｃ，２、ｑ_Ｃ，３、ｑ_Ｃ，４を以下のように最小二乗法を用いて概算する。ただし、Ａ_Ｃのエレメントの１つであるγ_Ｃは通常非常に小さいので、計算量と精度との兼ね合いを考慮して０として計算する。

ただし、

である。

なお、上記の式では、一つの平面上の観測点と撮像素子上の結像点の組が、一つの平面あたりｎ_１組あり、その平面をｎ_２回場所を変えてカメラで撮像した、観測点と結像点の組を全て使用している。

［Ｓｔｅｐ３．カメラ内部／外部パラメタ、レンズ歪み係数を最適化手法により高精度化］
次に、上記の式が最も小さくなるように、カメラ内部行列Ａ_Ｃ、レンズ歪み係数ｋ_Ｃ，１、ｋ_Ｃ，２、ｑ_Ｃ，１、ｑ_Ｃ，２、ｑ_Ｃ，３、ｑ_Ｃ，４、カメラ回転行列Ｒ_Ｃ，Ｂｊ、カメラ平行並進ベクトルｔ_Ｃ，Ｂｊを、たとえばLevenberg-Marquardt法などで最適化する。

ここで、ｎ_１は１枚の平面上にある観測点と撮像素子上の結像点の組の数、ｎ_２は平面の場所を変えて撮像した回数、Ｘ_{ｐＢｊ，ｉ}’は観測点の座標、ブレーヴェｍ_{ｐＣ，Ｂｊ，ｉ}は結像点の実測した座標、ｍ［．］は観測点Ｘ_{ｐＢｊ，ｉ}’をＡ_Ｃ、ｋ_Ｃ，１、ｋ_Ｃ，２、ｑ_Ｃ，１、ｑ_Ｃ，２、ｑ_Ｃ，３、ｑ_Ｃ，４、Ｒ_Ｃ，Ｂｊ、ｔ_Ｃ，Ｂｊを用いてカメラの撮像素子上に結像した結像点の座標を計算する投射関数である。この投射関数ｍ［．］の最後の引数Ｘ_{ｐＢｊ，ｉ}’は３次元空間中の点であり、ユーザが規定した世界座標系上のユーザが決めた点の座標であり、既知である。引数Ａ_Ｃ、ｋ_Ｃ，１、ｋ_Ｃ，２、ｑ_Ｃ，１、ｑ_Ｃ，２、ｑ_Ｃ，３、ｑ_Ｃ，４、Ｒ_Ｃ，Ｂｊ、ｔ_Ｃ，Ｂｊが、求めたいパラメタである。本実施の形態では投射関数ｍ［．］を用いてパラメタを求めているが、これを一般的にいうと、カメラ位置と３次元上の観測点で規定されるエピポーラ幾何学の拘束条件を用いてパラメタを求めるということになる。

投射関数ｍ［.］は、レンズ歪みがない場合に結像される座標ｘ_{ｐＣ，Ｂ１，１}＝［ｘ_{ｐＣ，Ｂ１，１}，ｙ_{ｐＣ，Ｂ１，１}］^Ｔを入力とし、レンズで歪んだ場合に結像される座標ブレーヴェｘ_{ｐＣ，Ｂ１，１}＝［ブレーヴェｘ_{ｐＣ，Ｂ１，１}，ブレーヴェｙ_{ｐＣ，Ｂ１，１}］^Ｔを出力するレンズ歪み関数Ｄ［ｋ_Ｃ，１、ｋ_Ｃ，２、ｑ_Ｃ，１、ｑ_Ｃ，２、ｑ_Ｃ，３、ｑ_Ｃ，４、ｘ_{ｐＣ，Ｂｊ，ｉ}］を用いて、

と表される。ここでブレーヴェｘ_{ｐＣ，Ｂｊ，ｉ}は、

によって計算され、レンズ歪みの関数Ｄ［ｋ_Ｃ，１、ｋ_Ｃ，２、ｑ_Ｃ，１、ｑ_Ｃ，２、ｑ_Ｃ，３、ｑ_Ｃ，４、ｘ_{ｐＣ，Ｂｊ，ｉ}］は（８．１０）式に従って計算される。

（８．１２）式の最適化は、上記実施例では、Levenberg-Marquardt法により行うとしたが、勾配法やNewton法などを用いてもよい。

また、実施例中では、ｎ_１は８以上の整数、ｎ_２は３以上の整数としているが、ｎ_１はＳｔｅｐ１−１で使用しているホモグラフィ行列Ｈ_Ｃ，Ｂ１を一意に決めるのに必要な観察点の組みの数、すなわちＬ_Ｃ，Ｂ１ ^ＴＬ_Ｃ，Ｂ１のランクが８以上（Ｌ_Ｃ，Ｂ１ ^ＴＬ_Ｃ，Ｂ１は９×９の行列）となるような観察点の組の数でなければならず、９組が望ましい。

また、ｎ_２はＳｔｅｐ１−２で求めているカメラ内部行列Ａ_Ｃの独立な変数５個を出すために必要なものであり、Ｖ_Ｃ ^ＴＶ_Ｃのランクが５以上（Ｖ_Ｃ ^ＴＶ_Ｃは６×６の行列）となるような撮影回数でなければならない。

以上のように求められた、最適化された歪み係数ｋ_Ｃ，１、ｋ_Ｃ，２、ｑ_Ｃ，１、ｑ_Ｃ，２、ｑ_Ｃ，３、ｑ_Ｃ，４を係数として持つ歪み関数Ｄ［.］や、最適化されたカメラパラメタ行列Ａ_Ｃおよび最適化された歪み係数ｋ_Ｃ，１、ｋ_Ｃ，２、ｑ_Ｃ，１、ｑ_Ｃ，２、ｑ_Ｃ，３、ｑ_Ｃ，４を係数として持つ投射関数ｍ［.］、および、それらの逆関数を用いることで、光学系に偏心に依存した歪みが存在する場合にも、画像の歪みの補正や、画像からの３次元復元の際の歪みの補正などの画像処理といったカメラ校正を行うことができる。

本実施の形態の効果としては、レンズの偏心に依存しない従来の歪みに加えて、レンズの偏心に依存する歪みも加えたレンズ歪みの補正を行うため、高精度なレンズ歪みの補正が可能という利点がある。

空間中の平面上の観測点とカメラ撮像素子上の結像点の関係を示す図である。 レンズの偏心とレンズ歪みの関係（偏心量の２つの表示方法）を示す図である。 レンズの偏心とレンズ歪みの関係（偏心量と歪みの関係）を示す図である。 レンズの偏心とレンズ歪みの関係（レンズ偏心の有無による歪み画像の違い）を示す図である。 処理のフローを示す図である。

符号の説明

１０１は空間中にある平面Ｂ１、１０２は平面Ｂ１上の観測点ｐ_Ｂ１，１、１１１はカメラＣ、１１２はレンズ、１１３はレンズ中心、１１４は撮像素子、１１５は撮像素子原点、１１６は観測点ｐ_Ｂ１，１の結像ｐ_{Ｃ，Ｂ１，１}である。

Claims (5)

1. カメラで撮影された画像から、カメラ光学系の歪みを補正するカメラ校正装置におけるカメラ校正方法であって、
   校正の際に使われるレンズ歪み関数が、偏心に依存するレンズ歪みをパラメタとして有し、
   撮像面上のレンズの光軸中心からの距離の３乗と撮像面の各軸方向への偏心量に比例したレンズ歪みと、撮像面上のレンズの光軸中心からの距離の５乗と撮像面の各軸方向への偏心量に比例したレンズ歪みとの和を含む式として表されること、
   を特徴とするカメラ校正方法。
2. カメラで撮影された画像から、カメラ光学系の歪みを補正するカメラ校正装置におけるカメラ校正方法であって、
   校正の際に使われるレンズ歪み関数が、偏心に依存するレンズ歪みをパラメタとして有し、
   一つの平面上にある観察点の座標と撮像面上の結像点の座標の組からレンズ歪みが無いとして、カメラ内部／外部パラメタを概算する第１のステップと、
   光学系の偏心に依存しない、撮像面上のレンズの光軸中心からの方向によらず、撮像面上のレンズの光軸中心からの距離のみにより歪み量が異なる歪みと、光学系の偏心に依存し、撮像面上のレンズの光軸中心からの方向と距離の両方により歪み量が異なる歪みの両方を足し合わせたレンズ歪み関数を用いて、３次元空間中の点のカメラ撮像面上の観測位置と、前記レンズ歪み関数を用いて計算したカメラ撮像上の位置との差の二乗が最も小さくなる最小二乗法で歪み係数を概算する第２のステップと、
   前記歪み関数と、カメラ位置と３次元上の観測点で規定されるエピポーラ幾何学の拘束条件を用いて、３次元空間中の点のカメラ撮像面上の観測位置と、前記レンズ歪み関数を用いて計算したカメラ撮像面上の位置との差が極小となるような最適化手法を用いてカメラ内部／外部パラメタとレンズ歪み係数を高精度化する第３のステップ、
   を持つことを特徴とするカメラ校正方法。
3. 請求項１または２に記載のカメラ校正方法をコンピュータに実行させるためのプログラム。
4. 請求項１または２に記載のカメラ校正方法をコンピュータに実行させるためのプログラムを格納した記憶媒体。
5. 請求項１または２に記載のカメラ校正方法を実行するカメラ校正装置。

Images