(57)【要約】【課題】偏心しているレンズを有するカメラにおいても、歪み補正を精度よく行うカメラ校正方法を提供する。【解決手段】カメラで撮影された画像から、カメラ光学系の歪みを補正するカメラ校正装置におけるカメラ校正方法であり、校正の際に使われるレンズ歪み関数が、偏心に依存するレンズ歪みをパラメタとして有する。Ｓｔｅｐ０で、観測点座標と結像点座標の組を入力する。Ｓｔｅｐ１で、レンズ歪みが無いとして、カメラ内部／外部パラメタを概算する。Ｓｔｅｐ２で、レンズ歪み係数を概算する。Ｓｔｅｐ３で、カメラ内部／外部パラメタ、レンズ歪み係数を最適化手法により高精度化する。この処理の流れは従来と同じであるが、本発明では、Ｓｔｅｐ２以降で求められるレンズ歪み係数の中に、レンズの偏心に依存した成分が新たに導入されている。【選択図】図３
2008262255.tif 000002
【特許請求の範囲】
【請求項１】
カメラで撮影された画像から、カメラ光学系の歪みを補正するカメラ校正装置におけるカメラ校正方法であって、
校正の際に使われるレンズ歪み関数が、偏心に依存するレンズ歪みをパラメタとして有すること、
を特徴とするカメラ校正方法。
【請求項２】
前記校正の際に使われるレンズ歪み関数が、
光学系の偏心に依存しない、撮像面上のレンズの光軸中心からの方向によらず、撮像面上のレンズの光軸中心からの距離のみより歪み量が異なる歪みと、光学系の偏心に依存し、撮像面上のレンズの光軸中心からの方向と距離の両方により歪み量が異なる歪みの和として表されること、
を特徴とする請求項１記載のカメラ校正方法。
【請求項３】
前記校正の際に使われるレンズ歪み関数が、
撮像面上のレンズの光軸中心からの距離の３乗と撮像面の各軸方向への偏心量に比例したレンズ歪みと、撮像面上のレンズの光軸中心からの距離の５乗と撮像面の各軸方向への偏心量に比例したレンズ歪みとの和として表されること、
を特徴とする請求項１もしくは２記載のカメラ校正方法。
【請求項４】
一つの平面上にある観察点の座標と撮像面上の結像点の座標の組からレンズ歪みが無いとして、カメラ内部／外部パラメタを概算する第一のステップと、
光学系の偏心に依存しない、撮像面上のレンズの光軸中心からの方向によらず、撮像面上のレンズの光軸中心からの距離のみにより歪み量が異なる歪みと、光学系の偏心に依存し、撮像面上のレンズの光軸中心からの方向と距離の両方により歪み量が異なる歪みの両方を足し合わせたレンズ歪み関数を用いて、
３次元空間中の点のカメラ撮像面上の観測位置と、前記レンズ歪み関数を用いて計算したカメラ撮像上の位置との差の二乗が最も小さくなる最小二乗法で歪み係数を概算する第２のステップと、
前記歪み関数と、カメラ位置と３次元上の観測点で規定されるエピポーラ幾何学の拘束条件を用いて、
３次元空間中の点のカメラ撮像面上の観測位置と、前記レンズ歪み関数を用いて計算したカメラ撮像面上の位置との差が極小となるような最適化手法を用いてカメラ内部／外部パラメタとレンズ歪み係数を高精度化する第３のステップ、
を持つことを特徴とする、請求項１ないし３のうちのいずれか1項に記載のカメラ校正方法。
【請求項５】
請求項１ないし４のうちいずれか１項に記載のカメラ校正方法をコンピュータに実行させるためのプログラム。
【請求項６】
請求項１ないし４のうちいずれか１項に記載のカメラ校正方法をコンピュータに実行させるためのプログラムを格納した記憶媒体。
【請求項７】
請求項１ないし４のうちいずれか１項に記載のカメラ校正方法を実行するカメラ校正装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、カメラで撮った複数枚の画像から、カメラの位置姿勢、カメラの内部パラメタ等を推定するカメラ校正技術に関する。
【背景技術】
【０００２】
カメラの光学系に歪みがない理想的な系を考えた場合、カメラ撮像素子上の点ｍ＝［ｕ，ｖ］Ｔと空間中の点Ｍ＝［Ｘ，Ｙ，Ｚ］Ｔとの間には、レンズの焦点距離、カメラ撮像素子の画像ピッチ(水平走査方向、垂直走査方向)、光学系の光軸と撮像素子との交点座標といったカメラ内部パラメタ、カメラの位置や姿勢といったカメラ外部パラメタを用いて、
【０００３】
【数１】
2008262255.tif 000003
【０００４】
ここで、αＣは、カメラ撮像素子の水平走査方向の画素ピッチと、カメラレンズの焦点距離の、積の逆数、βＣは、カメラ撮像素子の垂直走査方向の画素ピッチと、レンズの焦点距離の、積の逆数を、水平走査方向と垂直走査方向の成す角の正弦（ｓｉｎ）で割った値、γＣは、カメラ撮像素子の垂直方向の画素ピッチと、レンズ焦点距離の、積の逆数を、水平走査方向と垂直走査方向の成す角の正接（ｔａｎ）で割った値、ｕＣ，０は、光学系の光軸と撮像素子との交点のカメラ撮像素子の水平走査方向の座標、ｖＣ，０は、光学系の光軸と撮像素子との交点のカメラ撮像素子の垂直走査方向の座標、ＲＣはカメラ回転行列、ｔＣはカメラ平行並進ベクトル、のような関係が成り立つ。なお、上記の文章中に記載したｍ＝［ｕ，ｖ］Ｔ、Ｍ＝［Ｘ，Ｙ，Ｚ］Ｔの左辺の「ｍ」、「Ｍ」は本来は太字で記載すべきものであるが、文章中では太字が使えないため通常の文字で「ｍ」、「Ｍ」と記載した。以下の文章中の式で用いる文字についても同様である。
【０００５】
ところが実際には、光学系には一般に「樽型歪み」や「糸巻き型歪み」と呼ばれる歪曲歪みが存在しており、空間中の点Ｍの理想的なカメラ撮像素子上の射影点ｍと、実際にカメラ撮像素子上に結像される座標（すなわち、レンズで歪んだ場合に結像される座標）
【０００６】
【数２】
2008262255.tif 000004
【０００７】
とは一致しない。なお、数２のイメージデータに示した記号を文章中では「ブレーヴェｍ」と記載する。
【０００８】
また、一般にカメラ内部パラメタ、カメラ外部パラメタも厳密な値は既知ではないので、理想的なカメラ撮像素子上の射影点ｍから実際にカメラ撮像素子上に結像される座標ブレーヴェｍへの変換関数Ｄ［．］とともに、カメラ内部パラメタ、カメラ外部パラメタも、空間中の観測点の座標と撮像素子上の結像点の座標の組から推定する必要がある。
【０００９】
このような歪曲歪みのパラメタも含めたカメラパラメタの推定方法として、従来、以下のような方法が行われていた（たとえば、非特許文献１）。
【００１０】
図１は空間中にある平面上の観測点とカメラ撮像素子上の結像点の関係を示す図であり、１０１は空間中の平面Ｂ１、１０２は平面上の観測点ｐＢ１，１、１１１はカメラＣ、１１２はレンズ、１１３はレンズ中心、１１４は撮像素子、１１５は撮像素子原点、１１６は観測点ｐＢ１，１の結像点ｐＣ，Ｂ１，１、を示す。
【００１１】
カメラパラメタの推定は、たとえば、下記のステップで行われる。
Ｓｔｅｐ１．レンズ歪みが無いとして、カメラ内部／外部パラメタを概算。
Ｓｔｅｐ１−１．空間中の平面Ｂ１とカメラ撮像素子のホモグラフィ行列を概算。
Ｓｔｅｐ１−２．カメラ内部パラメタを概算。
Ｓｔｅｐ１−３．カメラ外部パラメタを概算。
Ｓｔｅｐ２．レンズ歪み係数を概算。
Ｓｔｅｐ３．カメラ内部／外部パラメタ、レンズ歪み係数を最適化手法により高精度化。
【００１２】
以下、詳細に説明する。
【００１３】
［Ｓｔｅｐ１．レンズ歪みが無いとして、カメラ内部／外部パラメタを概算］
［Ｓｔｅｐ１−１．空間中の平面Ｂ１とカメラ撮像素子のホモグラフィ行列を概算］
平面Ｂ１上の観測点ｐＢ１，１の座標を世界座標系で表し、ＸｐＢ１，１＝［ＸｐＢ１，１，ＹｐＢ１，１，ＺｐＢ１，１］Ｔのようにデカルト座標系で表記する。また、撮像素子上にある観測点ｐＢ１，１のカメラＣの撮像素子上の結像点ｐＣ，Ｂ１，１の座標を撮像素子のディジタル画像座標系で表し、ｍｐＣ，Ｂ１，１＝［ｕｐＣ，Ｂ１，１，ｖｐＣ，Ｂ１，１］Ｔのように表記する。
【００１４】
１０１の平面Ｂ１を世界座標系の座標［Ｘ，Ｙ，Ｚ］ＴのＺ＝０の平面とし、その座標を２次元座標ＸｐＢ１，１’＝［ＸｐＢ１，１，ＹｐＢ１，１］Ｔとしたとき、平面Ｂ１上の観測点ｐＢ１，１と結像点ｐＣ，Ｂ１，１のそれぞれの座標は、平面Ｂ１と撮像素子の平面で関連づけられる１つのホモグラフィ行列ＨＣ，Ｂ１＝［ｈＣ，Ｂ１，１  ｈＣ，Ｂ１，２  ｈＣ，Ｂ１，３］を使って次のように関係づけられる。ただし、ｈＣ，Ｂ１，１、ｈＣ，Ｂ１，２、ｈＣ，Ｂ１，３はそれぞれ３行１列の列ベクトルである。
【００１５】
【数３】
2008262255.tif 000005
【００１６】
ここで、ベクトルの上に付いている「〜」は同次座標であることを示し、チルダｍｐＣ，Ｂ１，１＝［ｕｐＣ，Ｂ１，１，ｖｐＣ，Ｂ１，１，１］Ｔ、チルダＸｐＢ１，１’＝［ＸｐＢ１，１，ＹｐＢ１，１，１］Ｔであり、ｓはスカラ値である。なお、数３のイメージに記載されているような「〜」が上についた「ｍ」等を文章中では「チルダｍ」等と記載する。
【００１７】
平面Ｂ１上の点の座標とそれに対応する撮像素子上の結像点の座標を（３．１）式に代入して方程式をつくり、ホモグラフィ行列ＨＣ，Ｂ１＝［ｈＣ，Ｂ１，１  ｈＣ，Ｂ１，２  ｈＣ，Ｂ１，３］を算出する。この方法としては、たとえば、下記の方法がある。まず、ｎ組の観測点の座標と結像点の座標の組を使って、次の行列ＬＣ，Ｂ１を作る。
【００１８】
【数４】
2008262255.tif 000006
【００１９】
ここで、０＝［０，０，０］Ｔである。
【００２０】
次にＬＣ，Ｂ１ＴＬＣ，Ｂ１の固有ベクトルを求め、その中の最も固有値が小さい固有ベクトルをｘＣ，Ｂ１＝［ｈＣ，Ｂ１，１Ｔ  ｈＣ，Ｂ１，２Ｔ  ｈＣ，Ｂ１，３Ｔ］とする。最後にこのベクトルｘＣ，Ｂ１の各エレメントを並び直すことによりホモグラフィ行列ＨＣ，Ｂ１＝［ｈＣ，Ｂ１，１  ｈＣ，Ｂ１，２  ｈＣ，Ｂ１，３］を得る。
【００２１】
［Ｓｔｅｐ１−２．カメラ内部パラメタを概算］
次に、ホモグラフィ行列ＨＣ，Ｂ１＝［ｈＣ，Ｂ１，１  ｈＣ，Ｂ１，２  ｈＣ，Ｂ１，３］のエレメントを次のように表記する。
【００２２】
【数５】
2008262255.tif 000007
【００２３】
ただし、ｉは１≦ｉ≦３の整数である。
【００２４】
上記表記を使って、次のようなベクトルを作る。
【００２５】
【数６】
2008262255.tif 000008
【００２６】
ただし、ｊは１≦ｊ≦３の整数である。
【００２７】
平面の位置をｎ２回変えてホモグラフィ行列ＨＣ，Ｂ１＝［ｈＣ，Ｂ１，１  ｈＣ，Ｂ１，２  ｈＣ，Ｂ１，３］を作り、それらから得たｎ２個のベクトルＶＣ，Ｂ１，ｉ，ｊ〜ＶＣ，Ｂｎ２，ｉ，ｊを使って、次の行列ＶＣを作る。
【００２８】
【数７】
2008262255.tif 000009
【００２９】
次にＶＣＴＶＣの固有ベクトルを求め、その中の最も固有値が小さい固有ベクトルをｂＣ＝［ＢＣ，１，１  ＢＣ，１，２  ＢＣ，２，２  ＢＣ，１，３  ＢＣ，２，３  ＢＣ，３，３］とする。
【００３０】
このｂＣを用いてカメラ内部行列ＡＣは以下のように概算される。
【００３１】
【数８】
2008262255.tif 000010
【００３２】
ただし、
【００３３】
【数９】
2008262255.tif 000011
【００３４】
ここで、αＣは、カメラ撮像素子の水平走査方向の画素ピッチと、レンズの焦点距離の、積の逆数、βＣは、カメラ撮像素子の垂直走査方向の画素ピッチと、レンズの焦点距離の、積の逆数を、水平走査方向と垂直走査方向の成す角の正弦（ｓｉｎ）で割ったもの、γＣは、カメラ撮像素子の水平走査方向の画素ピッチと、レンズの焦点距離の、積の逆数を、水平走査方向と垂直走査方向の成す角の正接（ｔａｎ）で割ったもの、ｕＣ，０は、光学系の光軸と撮像素子との交点のカメラ撮像素子の水平走査方向の座標、ｖＣ，０は、光学系の光軸と撮像素子との交点のカメラ撮像素子の垂直走査方向の座標、である。
【００３５】
［Ｓｔｅｐ１−３．カメラ外部パラメタを概算］
一方、カメラ外部行列ＧＣ，Ｂ１をカメラ外部パラメタをエレメントに持つカメラ回転行列ＲＣ，Ｂ１＝［ｒＣ，Ｂ１，１  ｒＣ，Ｂ１，２  ｒＣ，Ｂ１，３］とカメラ平行並進ベクトルｔＣ，Ｂ１からＧＣ，Ｂ１＝［ＲＣ，Ｂ１  ｔＣ，Ｂ１］と定義すると、カメラ外部行列ＧＣ，Ｂ１はこれまでの結果を用いて、
【００３６】
【数１０】
2008262255.tif 000012
【００３７】
と概算される。
【００３８】
［ｓｔｅｐ２．レンズ歪み係数を概算］
次に、レンズ歪みを反映するためのレンズ歪み係数を計算する。従来の方法においては、レンズで歪んだ場合に結像される座標ブレーヴェｘｐＣ，Ｂ１，１＝［ブレーヴェｘｐＣ，Ｂ１，１，ブレーヴェｙｐＣ，Ｂ１，１］Ｔは、レンズ歪みがない場合に結像される座標ｘｐＣ，Ｂ１，１＝［ｘｐＣ，Ｂ１，１，ｙｐＣ，Ｂ１，１］Ｔ、レンズ歪み係数ｋＣ，１、ｋＣ，２を用いて、
【００３９】
【数１１】
2008262255.tif 000013
【００４０】
としていた。ただし、本式はカメラＣのカメラ座標系で記述している。なお、数１０にイメージデータで示された（３．１０）式の左辺のｘ、ｙを文章中では「ブレーヴェｘ」、「ブレーヴェｙ」と記載する。
【００４１】
この変換を、レンズ歪みがない場合に結像される座標ｘｐＣ，Ｂ１，１＝［ｘｐＣ，Ｂ１，１，ｙｐＣ，Ｂ１，１］Ｔを入力とし、レンズで歪んだ場合に結像される座標ブレーヴェｘｐＣ，Ｂ１，１＝［ブレーヴェｘｐＣ，Ｂ１，１，ブレーヴェｙｐＣ，Ｂ１，１］Ｔを出力するレンズ歪み関数Ｄ［ｋＣ，１，ｋＣ，２，ｘｐＣ，Ｂ１，１］を導入し、ブレーヴェｘｐＣ，Ｂ１，１＝Ｄ［ｋＣ，１，ｋＣ，２，ｘｐＣ，Ｂ１，１］と表す。
【００４２】
上記の式をカメラ座標系ｘＣ＝［ｘＣ，ｙＣ］とディジタル座標系ｍＣ＝ＡＣｘＣの関係を使い、レンズ歪み係数ｋＣ，１、ｋＣ，２を以下のように算出する。ただし、ＡＣのエレメントの１つであるγＣは通常非常に小さいので、計算量と精度との兼ね合いを考慮して０として計算する。
【００４３】
【数１２】
2008262255.tif 000014
【００４４】
ただし、
【００４５】
【数１３】
2008262255.tif 000015
【００４６】
なお、上記の式では、一つの平面上の観測点と撮像素子上の結像点の組が、１つの平面あたりｎ１組あり、その平面をｎ２回場所を変えてカメラで撮像した、観測点と結像点の組を全て使用している。
【００４７】
［Ｓｔｅｐ３．カメラ内部／外部パラメタ、レンズ歪み係数を最適化手法により高精度化］
次に、下記の式が最も小さくなるように、カメラ内部行列ＡＣ、カメラ回転行列ＲＣ，Ｂｊ、カメラ平行並進ベクトルｔＣ，Ｂｊ、レンズ歪み係数ｋＣ，１、ｋＣ，２を最適化し、カメラ内部／外部パラメタおよびレンズ歪み係数を高精度化する。
【００４８】
【数１４】
2008262255.tif 000016
【００４９】
ここで、ｎ１は１枚の平面上にある観測点と撮像素子上の結像点の組の数、ｎ２は平面の場所を変えて撮像した回数、ＸｐＢｊ，ｉ’は観測点の座標、ブレーヴェｍｐＣ，Ｂｊ．ｉは結像点の実測した座標、ｍ［．］は観測点ＸｐＢｊ，ｉ’をＡＣ、ｋＣ，１、ｋＣ，２、ＲＣ，Ｂｊ、ｔＣ，Ｂｊを用いてカメラの撮像素子上に結像した結像点の座標を計算する投射関数であり、これまで求めてきた式を用いて、
【００５０】
【数１５】
2008262255.tif 000017
【００５１】
と表される。ここでチルダｘｐＣ，Ｂｊ，ｉは、
【００５２】
【数１６】
2008262255.tif 000018
【００５３】
によって計算される。（３．１２）式の最適化は、たとえば、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/.
【発明の開示】
【発明が解決しようとする課題】
【００５６】
しかしながら、従来の方法で求められたカメラ内部／外部パラメタおよび歪み係数を用いて実際に画像処理を行うと、レンズの歪みが残っている場合があった。特に近年多く使われるようになった、携帯電話のカメラなどの安価な光学系を用いたカメラで顕著であった。
【００５７】
上記に示すような従来のレンズ歪み補正方法は、放射方向のレンズ歪みはレンズの偏心がない場合の歪曲歪みに対応しているだけであり、レンズの光軸からの偏心がある場合の歪曲歪みについては考慮されておらず、偏心に依存した歪み量算出のためのカメラパラメタとして何を導入すればよいかの指針も明らかになっていなかったことに起因する。このため従来法では、一般のレンズが偏心しているカメラを使用した場合には歪み補正が十分なされない、という欠点があった。
【００５８】
この点について、図２を用いて説明する。レンズの偏心とは、光学系の基準軸に対してレンズが平行移動したり傾いたりすることにより、レンズ面の曲率中心が基準軸上からずれることを言う。図２（ａ）はこの様子を示したものであり、上図はずれ量を平行移動距離で表したものであり、下図は角度で表したものである。軸のずれがわずかな場合は、Ｅ＝ｌεであるので、平行移動Ｅと角度εは比例関係となるので、歪み係数を求める際には本質的な違いがなく、どちらで表しても良い。
【００５９】
図２（ｂ）は、偏心量をＥで表し、偏心した場合の歪み量について説明する図である。偏心した場合の歪み量は、無歪み時の結像位置ｌに加えて、偏心に依存しない場合の放射方向歪み（Ｖ３‖ｌ‖２＋Ｖ５‖ｌ‖４）と、偏心に依存する放射方向歪み（Ｖ３Ｅ＋Ｖ５Ｅ‖ｌ‖２）（Ｅ・ｌ）ｌがある。従来は、前者の歪みのみを考慮していたが、後者の歪みもあるため、従来の方法では偏心がある場合は後者の歪みの補正が十分なされない。この理由としては、後者は（Ｅ・ｌ）の項を持つため基準軸からの方向によってその歪みの大きさが変わるという性質を持つが、従来考慮していた前者の歪みは基準軸からの方向により歪みの大きさが変わらないため、前者の歪みを考慮するだけでは、後者の歪みは十分補正できない。
【００６０】
基準軸からの方向により歪む量が異なることについて、結像する像の形を用いて説明する。図２（ｃ）にレンズの偏心がある場合と無い場合とでの歪み画像の違いを示す。点線は無歪み時の画像であり、一点鎖線はレンズの偏心が無い場合の歪み画像であり、実線はレンズの偏心がある場合の歪み画像である。偏心がなければ、光学系の基準軸を中心にどの角度にも同じ歪みが出て、一点鎖線のような画像となる。偏心があった場合には、矢印に示す偏心の方向に像が伸縮する。従来は、一点鎖線で示した偏心が無い場合の方向性のない歪みのみを補正対象としているが、レンズの偏心がある一般的な場合では、レンズの偏心による実線で示したレンズの偏心方向に依存した歪みが出ており、この歪みは従来方法では十分に補正することはできない。
【００６１】
本発明は、上記に鑑みてなされたものであり、その目的とするところは、偏心しているレンズを有するカメラにおいても、歪み補正を精度よく行うカメラ校正方法を提供することにある。
【課題を解決するための手段】
【００６２】
本願において開示される発明のうち、代表的なものの概要を簡単に説明すれば、以下の通りである。
【００６３】
第１の発明は、カメラで撮影された画像から、カメラ光学系の歪みを補正するカメラ校正装置におけるカメラ校正方法であって、校正の際に使われるレンズ歪み関数が、偏心に依存するレンズ歪みをパラメタとして有すること、を特徴とする。
【００６４】
第２の発明は、第１の発明において、前記校正の際に使われるレンズ歪み関数が、光学系の偏心に依存しない、撮像面上のレンズの光軸中心からの方向によらず、撮像面上のレンズの光軸中心からの距離のみより歪み量が異なる歪みと、光学系の偏心に依存し、撮像面上のレンズの光軸中心からの方向と距離の両方により歪み量が異なる歪みの和として表されること、を特徴とする。
【００６５】
第３の発明は、第１または第２の発明において、前記校正の際に使われるレンズ歪み関数が、撮像面上のレンズの光軸中心からの距離の３乗と撮像面の各軸方向への偏心量に比例したレンズ歪みと、撮像面上のレンズの光軸中心からの距離の５乗と撮像面の各軸方向への偏心量に比例したレンズ歪みとの和として表されること、を特徴とする。
【００６６】
第４の発明は、第１〜３に記載の発明において、一つの平面上にある観察点の座標と撮像面上の結像点の座標の組からレンズ歪みが無いとして、カメラ内部／外部パラメタを概算する第一のステップと、光学系の偏心に依存しない、撮像面上のレンズの光軸中心からの方向によらず、撮像面上のレンズの光軸中心からの距離のみにより歪み量が異なる歪みと、光学系の偏心に依存し、撮像面上のレンズの光軸中心からの方向と距離の両方により歪み量が異なる歪みの両方を足し合わせたレンズ歪み関数を用いて、３次元空間中の点のカメラ撮像面上の観測位置と、前記レンズ歪み関数を用いて計算したカメラ撮像上の位置との差の二乗が最も小さくなる最小二乗法で歪み係数を概算する第２のステップと、前記歪み関数と、カメラ位置と３次元上の観測点で規定されるエピポーラ幾何学の拘束条件を用いて、３次元空間中の点のカメラ撮像面上の観測位置と、前記レンズ歪み関数を用いて計算したカメラ撮像面上の位置との差が極小となるような最適化手法を用いてカメラ内部／外部パラメタとレンズ歪み係数を高精度化する第３のステップ、を持つことを特徴とする。
【００６７】
第５の発明は、第１〜４の発明のカメラ校正方法をコンピュータに実行させるためのプログラムである。
【００６８】
第６の発明は、第１〜４の発明のカメラ校正方法をコンピュータに実行させるためのプログラムを格納した記憶媒体である。
【００６９】
第７の発明は、請求項１〜４の発明のカメラ校正方法を実行するカメラ校正装置である。
【００７０】
以下、本発明の理解を助けるために、より具体的に説明するが、これにより本発明が限定されるものではない。
【発明の効果】
【００７１】
本発明により、レンズの偏心に依存しない従来の歪みに加えて、レンズの偏心に依存する歪みも加えたレンズ歪みの補正を行うため、高精度なレンズ歪みの補正が可能となる。
【発明を実施するための最良の形態】
【００７２】
以下、図面を用いて本発明の実施の形態を説明する。
【００７３】
図１は空間中にある平面上の観測点とカメラ撮像素子上の結像点の関係を示す図であり、１０１は空間中の平面Ｂ１、１０２は平面上の観測点ｐＢ１，１、１１１はカメラＣ、１１２はレンズ、１１３はレンズ中心、１１４は撮像素子、１１５は撮像素子原点、１１６は観測点ｐＢ１，１の結像点ｐＣ，Ｂ１，１、を示す。
【００７４】
図３に本実施形態のカメラ校正方法の処理のフローを示す。なお、本明細書においては、カメラ内部／外部パラメタおよびレンズ歪み係数等のパラメタの算出方法をカメラ校正方法と呼ぶ。
【００７５】
カメラ校正に用いられるカメラ内部／外部パラメタおよびレンズ歪み係数は、以下のように求められる。以下の処理はカメラ校正装置（図示していない）によって実行することができる。カメラ校正装置は記憶装置に記憶されたプログラムとコンピュータによって構成することができ、また、そのプログラムの一部または全部に代えてハードウェアで構成することもできる。
Ｓｔｅｐ０．観測点座標と結像点座標の組の入力。
Ｓｔｅｐ１．レンズ歪みが無いとして、カメラ内部／外部パラメタを概算。
Ｓｔｅｐ１−１．空間中の平面Ｂ１とカメラ撮像素子のホモグラフィ行列を概算。
Ｓｔｅｐ１−２．カメラ内部パラメタを概算。
Ｓｔｅｐ１−３．カメラ外部パラメタを概算。
Ｓｔｅｐ２．レンズ歪み係数を概算。
Ｓｔｅｐ３．カメラ内部／外部パラメタ、レンズ歪み係数を最適化手法により高精度化。
【００７６】
処理の流れは従来と同じであるが、Ｓｔｅｐ２以降で求められるレンズ歪み係数の中に、レンズの偏心に依存した成分が新たに導入されている。
【００７７】
カメラ内部／外部パラメタおよびレンズ歪み係数（歪曲歪みの係数およびレンズ偏心の歪み係数）を求めるために、一平面上にあるｎ１個（ｎ１は８以上の整数）の観察点を、カメラと平面の相対位置を変えてｎ２回（ｎ２は３以上の整数）撮影したときの、空間中にある観察点の座標とカメラ撮像素子上の結像点の座標の組を入力として用いる。
【００７８】
以下、カメラ内部／外部パラメタおよびレンズ歪み係数の求め方を詳細に説明する。
［Ｓｔｅｐ０．観測点座標と結像点座標の組の入力］
平面Ｂ１上の観測点の座標と、それが撮像素子上で結像している点の座標の組を入力する。このとき、平面Ｂ１上の観測点の座標は、世界座標系で表し、また、撮像素子上で結像している点の座標は、カメラ座標軸で表す。
【００７９】
［Ｓｔｅｐ１．レンズ歪みが無いとして、カメラ内部／外部パラメタを概算］
［Ｓｔｅｐ１−１．空間中の平面Ｂ１とカメラ撮像素子のホモグラフィ行列を概算］
平面Ｂ１上の観測点ｐＢ１，１の座標を世界座標系で表し、ＸｐＢ１，１＝［ＸｐＢ１，１，ＹｐＢ１，１，ＺｐＢ１，１］Ｔのようにデカルト座標系で表記する。また、撮像素子上にある観測点ｐＢ１，１のカメラＣの撮像素子上の結像点ｐＣ，Ｂ１，１の座標を撮像素子のディジタル画像座標系で表し、ｍｐＣ，Ｂ１，１＝［ｕｐＣ，Ｂ１，１，ｖｐＣ，Ｂ１，１］Ｔのように表記する。
【００８０】
１０１の平面Ｂ１を世界座標系の座標［Ｘ，Ｙ，Ｚ］ＴのＺ＝０の平面とし、その座標を２次元座標ＸｐＢ１，１’＝［ＸｐＢ１，１，ＹｐＢ１，１］Ｔとしたとき、平面Ｂ１上の観測点ｐＢ１，１と結像点ｐＣ，Ｂ１，１のそれぞれの座標は、平面Ｂ１と撮像素子の平面で関連づけられる１つのホモグラフィ行列ＨＣ，Ｂ１＝［ｈＣ，Ｂ１，１  ｈＣ，Ｂ１，２  ｈＣ，Ｂ１，３］を使って次のように関係づけられる。ただし、ｈＣ，Ｂ１，１、ｈＣ，Ｂ１，２、ｈＣ，Ｂ１，３はそれぞれ３行１列の列ベクトルである。
【００８１】
【数１７】
2008262255.tif 000019
【００８２】
ここで、ベクトルの上に付いている「〜」は同次座標であることを示し、チルダｍｐＣ，Ｂ１，１＝［ｕｐＣ，Ｂ１，１，ｖｐＣ，Ｂ１，１，１］Ｔ、チルダＸｐＢ１，１’＝［ＸｐＢ１，１，ＹｐＢ１，１，１］Ｔであり、ｓはスカラ値である。
【００８３】
平面Ｂ１の点の座標とそれに対応する撮像素子上の結像点の座標を（８．１）式に代入して方程式をつくり、ホモグラフィ行列ＨＣ，Ｂ１＝［ｈＣ，Ｂ１，１  ｈＣ，Ｂ１，２  ｈＣ，Ｂ１，３］を算出する。この方法としては、たとえば、下記の方法がある。まず、ｎ１組（ただしｎ１は８以上の整数）の観測点の座標と結像点の座標の組を使って、次の行列ＬＣ，Ｂ１を作る。
【００８４】
【数１８】
2008262255.tif 000020
【００８５】
ここで、０＝［０，０，０］Ｔである。
【００８６】
次にＬＣ，Ｂ１ＴＬＣ，Ｂ１の固有ベクトルを求め、その中の最も固有値が小さい固有ベクトルをｘＣ，Ｂ１＝［ｈＣ，Ｂ１，１Ｔ  ｈＣ，Ｂ１，２Ｔ  ｈＣ，Ｂ１，３Ｔ］とする。最後にこのベクトルｘＣ，Ｂ１の各エレメントを並び直すことによりホモグラフィ行列ＨＣ，Ｂ１＝［ｈＣ，Ｂ１，１  ｈＣ，Ｂ１，２  ｈＣ，Ｂ１，３］を得る。
【００８７】
［Ｓｔｅｐ１−２．カメラ内部パラメタを概算］
次に、ホモグラフィ行列ＨＣ，Ｂ１＝［ｈＣ，Ｂ１，１  ｈＣ，Ｂ１，２  ｈＣ，Ｂ１，３］のエレメントを次のように表記する。
【００８８】
【数１９】
2008262255.tif 000021
【００８９】
ただしｉは１≦ｉ≦３の整数である。
【００９０】
上記表記を使って、次のようなベクトルを作る。
【００９１】
【数２０】
2008262255.tif 000022
【００９２】
ただしｊは１≦ｊ≦３の整数である。
【００９３】
平面の位置をｎ２回（ただしｎ２は３以上の整数）変えてホモグラフィ行列ＨＣ，Ｂ１＝［ｈＣ，Ｂ１，１  ｈＣ，Ｂ１，２  ｈＣ，Ｂ１，３］を作り、それらから得たｎ２個のベクトルｖＣ，Ｂ１，ｉ，ｊ〜ｖＣ，Ｂｍ，ｉ，ｊを使って、次の行列ＶＣを作る。
【００９４】
【数２１】
2008262255.tif 000023
【００９５】
次にＶＣＴＶＣの固有ベクトルを求め、その中の最も固有値が小さい固有ベクトルをｂＣ＝［ＢＣ，１，１  ＢＣ，１，2  ＢＣ，2，2  ＢＣ，１，3  ＢＣ，2，3  ＢＣ，3，3］とする。
【００９６】
このｂＣを用いてカメラ内部行列ＡＣは以下のように概算される。
【００９７】
【数２２】
2008262255.tif 000024
【００９８】
ただし、
【００９９】
【数２３】
2008262255.tif 000025
【０１００】
ここで、αＣは、カメラ撮像素子の水平走査方向の画素ピッチと、レンズの焦点距離の、積の逆数、βＣは、カメラ撮像素子の垂直走査方向の画素ピッチと、レンズの焦点距離の、積の逆数を、水平走査方向と垂直走査方向の成す角の正弦（ｓｉｎ）で割ったもの、γＣは、カメラ撮像素子の水平走査方向の画素ピッチと、レンズの焦点距離の、積の逆数を、水平走査方向と垂直走査方向の成す角の正接（ｔａｎ）で割ったもの、ｕＣ，０は、光学系の光軸と撮像素子との交点のカメラ撮像素子の水平走査方向の座標、ｖＣ，０は、光学系の光軸と撮像素子との交点のカメラ撮像素子の垂直走査方向の座標、である。
【０１０１】
［Ｓｔｅｐ１−３．カメラ外部パラメタを概算］
一方、カメラ外部行列ＧＣ，Ｂ１を、カメラ外部パラメタをエレメントに持つカメラ回転行列ＲＣ，Ｂ１とカメラ平行並進ベクトルｔＣ，Ｂ１からＧＣ，Ｂ１＝［ＲＣ，Ｂ１ｔＣ，Ｂ１］と定義すると、カメラ外部行列ＧＣ，Ｂ１はこれまでの結果を用いて
【０１０２】
【数２４】
2008262255.tif 000026
【０１０３】
ただし、
【０１０４】
【数２５】
2008262255.tif 000027
【０１０５】
と概算される。
【０１０６】
［Ｓｔｅｐ２．レンズ歪み係数を概算］
次に、レンズ歪みを反映するためのレンズ歪み係数を計算する。ここで考慮する歪みは、レンズ歪みの中で一般に歪み量の大きい、放射方向歪み（歪曲歪み）、および、レンズの曲率中心がレンズ基準軸からずれている状態、つまり、偏心している状態における歪みである。
【０１０７】
図２にレンズの偏心と歪みの関係を示す。偏心には図２（ａ）に示すように、２種類の表示方法があり、１つは図２（ａ）の上図に示すように並進移動Ｅによる表記、もう１つは図２（ａ）の下図に示すように角度εによる表記である。ずれ量が小さい場合は、曲率半径Ｒを用いてＥ＝Ｒεとなり、平行移動Ｅと角度εは比例関係となるので、歪み係数を求める際にはどちらを用いても本質的には違いがなく、どちらを用いても良い。一般に、レンズの偏心による歪みは、平行移動Ｅと角度εを比例関係として扱ってよい程度である。
【０１０８】
図２（ｂ）はＥを用いて歪みを示した図である。無歪み時の結像位置ｌに対して、
偏心に依存しない歪み（Ｖ３‖ｌ‖２＋Ｖ５‖ｌ‖４）と、それに加えて偏心に依存する歪み（Ｖ３Ｅ＋Ｖ５Ｅ‖ｌ‖２）（Ｅ・ｌ）ｌがあり、それらを加えた
【０１０９】
【数２６】
2008262255.tif 000028
【０１１０】
が全歪みとなる。
【０１１１】
これらをカメラＣのカメラ座標系の各エレメントで書き下したものが下記式（８.１０）である。ただし、係数Ｖ３、Ｖ５、Ｖ３Ｅ、Ｖ５Ｅ、Ｅ＝［Ｅｘ  Ｅｙ］はそのまま使わず、後述する（８．１１）式によって行列演算しやすいように、各項に乗じる係数を適宜組み合わせ、それらをｋＣ，１、ｋＣ，２、ｑＣ，１、ｑＣ，２、ｑＣ，３、ｑＣ，４という係数に置き換えた。
【０１１２】
【数２７】
2008262255.tif 000029
【０１１３】
ただし、
【０１１４】
【数２８】
2008262255.tif 000030
【０１１５】
である。
【０１１６】
ここで、ブレーヴェｘｐＣ，Ｂ１，１＝［ブレーヴェｘｐＣ，Ｂ１，１，ブレーヴェｙｐＣ，Ｂ１，１］Ｔはレンズで歪んだ場合に結像される座標、ｘｐＣ，Ｂ１，１＝［ｘｐＣ，Ｂ１，１，ｙｐＣ，Ｂ１，１］Ｔはレンズ歪みがない場合に結像される座標、ｋＣ，１はレンズの偏心に依存しない放射方向のレンズ歪み係数で、レンズ中心からの距離の３乗に比例した歪みの係数、ｋＣ，２はレンズの偏心に依存しない放射方向のレンズ歪み係数で、レンズ中心からの距離の５乗に比例した歪みの係数、ｑＣ，１はレンズの偏心に依存した放射方向のレンズ歪み係数で、レンズ中心からの距離の３乗と、撮像素子の水平走査方向の偏心量に依存した歪みの係数、ｑＣ，２はレンズの偏心に依存した放射方向のレンズ歪み係数で、レンズ中心からの距離の５乗と、撮像素子の水平走査方向の偏心量に依存した歪みの係数、ｑＣ，３はレンズの偏心に依存した放射方向のレンズ歪み係数で、レンズ中心からの距離の３乗と、撮像素子の垂直走査方向の偏心量に依存した歪みの係数、ｑＣ，４はレンズの偏心に依存した放射方向のレンズ歪み係数で、レンズ中心からの距離の５乗と、撮像素子の垂直走査方向の偏心量に依存した歪みの係数、である。
【０１１７】
（８−１０）式の上の２つの式のそれぞれの第２項目は偏心に関係なく現れる歪み量を表しており、第３項目、第４項目は偏心に依存する歪み量を表している。
【０１１８】
上記の式をカメラ座標系ｘＣ＝［ｘＣ，ｙＣ］とディジタル座標系ｍＣ＝ＡＣｘＣの関係を使い、レンズ歪み係数ｋＣ，１、ｋＣ，２、ｑＣ，１、ｑＣ，２、ｑＣ，３、ｑＣ，４を以下のように最小二乗法を用いて概算する。ただし、ＡＣのエレメントの１つであるγＣは通常非常に小さいので、計算量と精度との兼ね合いを考慮して０として計算する。
【０１１９】
【数２９】
2008262255.tif 000031
【０１２０】
ただし、
【０１２１】
【数３０】
2008262255.tif 000032
【０１２２】
である。
【０１２３】
なお、上記の式では、一つの平面上の観測点と撮像素子上の結像点の組が、一つの平面あたりｎ１組あり、その平面をｎ２回場所を変えてカメラで撮像した、観測点と結像点の組を全て使用している。
【０１２４】
［Ｓｔｅｐ３．カメラ内部／外部パラメタ、レンズ歪み係数を最適化手法により高精度化］
次に、上記の式が最も小さくなるように、カメラ内部行列ＡＣ、レンズ歪み係数ｋＣ，１、ｋＣ，２、ｑＣ，１、ｑＣ，２、ｑＣ，３、ｑＣ，４、カメラ回転行列ＲＣ，Ｂｊ、カメラ平行並進ベクトルｔＣ，Ｂｊを、たとえばLevenberg-Marquardt法などで最適化する。
【０１２５】
【数３１】
2008262255.tif 000033
【０１２６】
ここで、ｎ１は１枚の平面上にある観測点と撮像素子上の結像点の組の数、ｎ２は平面の場所を変えて撮像した回数、ＸｐＢｊ，ｉ’は観測点の座標、ブレーヴェｍｐＣ，Ｂｊ，ｉは結像点の実測した座標、ｍ［．］は観測点ＸｐＢｊ，ｉ’をＡＣ、ｋＣ，１、ｋＣ，２、ｑＣ，１、ｑＣ，２、ｑＣ，３、ｑＣ，４、ＲＣ，Ｂｊ、ｔＣ，Ｂｊを用いてカメラの撮像素子上に結像した結像点の座標を計算する投射関数である。この投射関数ｍ［．］の最後の引数ＸｐＢｊ，ｉ’は３次元空間中の点であり、ユーザが規定した世界座標系上のユーザが決めた点の座標であり、既知である。引数ＡＣ、ｋＣ，１、ｋＣ，２、ｑＣ，１、ｑＣ，２、ｑＣ，３、ｑＣ，４、ＲＣ，Ｂｊ、ｔＣ，Ｂｊが、求めたいパラメタである。本実施の形態では投射関数ｍ［．］を用いてパラメタを求めているが、これを一般的にいうと、カメラ位置と３次元上の観測点で規定されるエピポーラ幾何学の拘束条件を用いてパラメタを求めるということになる。
【０１２７】
投射関数ｍ［.］は、レンズ歪みがない場合に結像される座標ｘｐＣ，Ｂ１，１＝［ｘｐＣ，Ｂ１，１，ｙｐＣ，Ｂ１，１］Ｔを入力とし、レンズで歪んだ場合に結像される座標ブレーヴェｘｐＣ，Ｂ１，１＝［ブレーヴェｘｐＣ，Ｂ１，１，ブレーヴェｙｐＣ，Ｂ１，１］Ｔを出力するレンズ歪み関数Ｄ［ｋＣ，１、ｋＣ，２、ｑＣ，１、ｑＣ，２、ｑＣ，３、ｑＣ，４、ｘｐＣ，Ｂｊ，ｉ］を用いて、
【０１２８】
【数３２】
2008262255.tif 000034
【０１２９】
と表される。ここでブレーヴェｘｐＣ，Ｂｊ，ｉは、
【０１３０】
【数３３】
2008262255.tif 000035
【０１３１】
によって計算され、レンズ歪みの関数Ｄ［ｋＣ，１、ｋＣ，２、ｑＣ，１、ｑＣ，２、ｑＣ，３、ｑＣ，４、ｘｐＣ，Ｂｊ，ｉ］は（８．１０）式に従って計算される。
【０１３２】
（８．１２）式の最適化は、上記実施例では、Levenberg-Marquardt法により行うとしたが、勾配法やNewton法などを用いてもよい。
【０１３３】
また、実施例中では、ｎ１は８以上の整数、ｎ２は３以上の整数としているが、ｎ１はＳｔｅｐ１−１で使用しているホモグラフィ行列ＨＣ，Ｂ１を一意に決めるのに必要な観察点の組みの数、すなわちＬＣ，Ｂ１ＴＬＣ，Ｂ１のランクが８以上（ＬＣ，Ｂ１ＴＬＣ，Ｂ１は９×９の行列）となるような観察点の組の数でなければならず、９組が望ましい。
【０１３４】
また、ｎ２はＳｔｅｐ１−２で求めているカメラ内部行列ＡＣの独立な変数５個を出すために必要なものであり、ＶＣＴＶＣのランクが５以上（ＶＣＴＶＣは６×６の行列）となるような撮影回数でなければならない。
【０１３５】
以上のように求められた、最適化された歪み係数ｋＣ，１、ｋＣ，２、ｑＣ，１、ｑＣ，２、ｑＣ，３、ｑＣ，４を係数として持つ歪み関数Ｄ［.］や、最適化されたカメラパラメタ行列ＡＣおよび最適化された歪み係数ｋＣ，１、ｋＣ，２、ｑＣ，１、ｑＣ，２、ｑＣ，３、ｑＣ，４を係数として持つ投射関数ｍ［.］、および、それらの逆関数を用いることで、光学系に偏心に依存した歪みが存在する場合にも、画像の歪みの補正や、画像からの３次元復元の際の歪みの補正などの画像処理といったカメラ校正を行うことができる。
【０１３６】
本実施の形態の効果としては、レンズの偏心に依存しない従来の歪みに加えて、レンズの偏心に依存する歪みも加えたレンズ歪みの補正を行うため、高精度なレンズ歪みの補正が可能という利点がある。
【図面の簡単な説明】
【０１３７】
【図１】空間中の平面上の観測点とカメラ撮像素子上の結像点の関係を示す図である。
【図２（ａ）】レンズの偏心とレンズ歪みの関係（偏心量の２つの表示方法）を示す図である。
【図２（ｂ）】レンズの偏心とレンズ歪みの関係（偏心量と歪みの関係）を示す図である。
【図２（ｃ）】レンズの偏心とレンズ歪みの関係（レンズ偏心の有無による歪み画像の違い）を示す図である。
【図３】処理のフローを示す図である。
【符号の説明】
【０１３８】
１０１は空間中にある平面Ｂ１、１０２は平面Ｂ１上の観測点ｐＢ１，１、１１１はカメラＣ、１１２はレンズ、１１３はレンズ中心、１１４は撮像素子、１１５は撮像素子原点、１１６は観測点ｐＢ１，１の結像ｐＣ，Ｂ１，１である。
【図１】
2008262255.tif 000036
【図２（ａ）】
2008262255.tif 000037
【図２（ｂ）】
2008262255.tif 000038
【図２（ｃ）】
2008262255.tif 000039
【図３】
2008262255.tif 000040