(57)【要約】【課題】レンズ歪み係数を含む内部パラメータと外部パラメータとを高精度に推定し得る、カメラパラメータ推定装置、カメラパラメータ推定方法、およびプログラムを提供する。【解決手段】カメラパラメータ推定装置１０は、カメラによる対象物の撮影時におけるカメラパラメータを推定するための装置である。カメラパラメータ推定装置１０は、対象物に関する５組以上の３次元点と、３次元点それぞれに対応する画像上の２次元点との変換を表す誤差関数を、カメラパラメータの依存関係に基づいて、複数の部分問題に分割し、複数の部分問題それぞれにおいて停留点を算出する、パラメータ算出部２０と、誤差関数によって得られる誤差を最小とする停留点を求め、求めた停留点を最適なカメラパラメータとして出力する、パラメータ最適化部３０とを備えている。【選択図】図１
2017163386.tif 000002
【特許請求の範囲】
【請求項１】
カメラによる対象物の撮影時におけるカメラパラメータを推定するための装置であって、
前記対象物に関する５組以上の３次元点と、前記３次元点それぞれに対応する画像上の２次元点との変換を表す誤差関数を、前記カメラパラメータの依存関係に基づいて、複数の部分問題に分割し、前記複数の部分問題それぞれにおいて停留点を算出する、パラメータ算出部と、
前記誤差関数によって得られる誤差を最小とする停留点を求め、求めた前記停留点を最適なカメラパラメータとして出力する、パラメータ最適化部と、
を備えている、ことを特徴とするカメラパラメータ推定装置。
【請求項２】
前記カメラパラメータが、前記カメラの外部に関する外部パラメータを含み、
前記パラメータ算出部が、
前記複数の部分問題のうち、前記外部パラメータの一部のみに依存する１つの部分問題について、前記停留点を算出し、
更に、算出した前記停留点と、前記３次元点及び前記２次元点とに基づいて、前記停留点が未だ算出されていない前記部分問題について、前記停留点を算出する、
請求項１に記載のカメラパラメータ推定装置。
【請求項３】
前記パラメータ最適化部が、前記パラメータ算出部によって算出された前記停留点と、前記３次元点及び前記２次元点とに基づいて、前記誤差関数を最適化することによって、前記誤差を最小とする停留点を求める、
請求項１または２に記載のカメラパラメータ推定装置。
【請求項４】
前記外部パラメータが、前記カメラの３次元区間における位置及び姿勢のいずれか一方または両方を含み、
前記カメラパラメータが、前記外部パラメータに加えて、前記カメラの内部に関する内部パラメータを含み、前記内部パラメータは、前記カメラにおける、レンズの焦点距離、光学中心、アスペクト比、せん断係数、及びレンズ歪み係数のうち、少なくとも１つを含む、
請求項２に記載のカメラパラメータ推定装置。
【請求項５】
カメラによる対象物の撮影時におけるカメラパラメータを推定するための方法であって、
（ａ）前記対象物に関する５組以上の３次元点と、前記３次元点それぞれに対応する画像上の２次元点との変換を表す誤差関数を、前記カメラパラメータの依存関係に基づいて、複数の部分問題に分割し、前記複数の部分問題それぞれにおいて停留点を算出する、ステップと、
（ｂ）前記誤差関数によって得られる誤差を最小とする停留点を求め、求めた前記停留点を最適なカメラパラメータとして出力する、ステップと、
を有する、ことを特徴とするカメラパラメータ推定方法。
【請求項６】
コンピュータによって、カメラによる対象物の撮影時におけるカメラパラメータを推定するためのプログラムであって、
前記コンピュータに、
（ａ）前記対象物に関する５組以上の３次元点と、前記３次元点それぞれに対応する画像上の２次元点との変換を表す誤差関数を、前記カメラパラメータの依存関係に基づいて、複数の部分問題に分割し、前記複数の部分問題それぞれにおいて停留点を算出する、ステップと、
（ｂ）前記誤差関数によって得られる誤差を最小とする停留点を求め、求めた前記停留点を最適なカメラパラメータとして出力する、ステップと、
を実行させるプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、カメラの外部パラメータと内部パラメータとを推定する、カメラパラメータ推定装置及びカメラパラメータ推定方法に関し、更に、これらを実現するためのプログラムに関する。
【背景技術】
【０００２】
既知の３次元点を観測した画像からカメラパラメータを推定する技術は、ロボットの自己位置推定、および画像合成において、重要な要素技術である（例えば、非特許文献１〜５参照）。また、カメラパラメータには、外部パラメータと、内部パラメータとの２種類が存在する。外部パラメータは、カメラの３次元空間における位置と姿勢とを表すパラメータであり、内部パラメータは、レンズの焦点距離、光学中心、アスペクト比、せん断係数、レンズ歪み係数などのパラメータである。以下では、既知の３次元点と、この３次元点を画像上で観測することによって得られた２次元点とのペアを、単に「対応点」と表記する。
【０００３】
非特許文献１と非特許文献２とは、共に、表面模様の３次元点が既知の平面、または立体のパターンを複数の角度から撮影して得られた画像を用いて、外部パラメータと内部パラメータとを計算する方法を開示している。これらの非特許文献に開示された方法では、まずレンズ歪み係数がすべてゼロに仮定され、姿勢に関する拘束条件を無視した上で内部パラメータと外部パラメータとが計算される。そして、計算された内部パラメータおよび外部パラメータを初期値として、ニュートン法などの非線形最適化が行われる。
【０００４】
また、非特許文献３は、非特許文献１または非特許文献２に開示された方法を利用して、予め内部パラメータを計算しておき、３組以上の対応点から外部パラメータのみを計算する方法を開示している。具体的には、非特許文献３に開示された方法は、３次元点と２次元点との射影関係に基づく誤差関数を多項式で表現し、誤差関数の停留点である勾配のゼロ点をすべて計算することで、局所解を回避して大域的に最適な唯一の解を推定する。
【０００５】
更に、非特許文献４は、内部パラメータのうち焦点距離のみを未知として、４組以上の対応点を用いて、非特許文献３に開示された方法と同様にして、多項式表現による大域的最適解を直接計算する方法を開示している。
【０００６】
また、非特許文献５と非特許文献６とは、それぞれ、４組または５組の対応点を入力として、内部パラメータのうち焦点距離とレンズ歪み係数を未知として、これらの内部パラメータを外部パラメータとを共に算出する方法を開示している。
【先行技術文献】
【非特許文献】
【０００７】
【非特許文献１】Z. Zhang, "Flexible camera calibration by viewing a plane from unknown orientations," ICCV 1999.
【非特許文献２】R. Tsai, "A versatile camera calibration technique for high-accuracy 3D machine vision metrology using off-the-shelf TV cameras and lenses," Robotics and Automation, IEEE Journal of 3.4 (1987): 323-344.
【非特許文献３】G. Nakano, “Globally Optimal DLS Method for PnP Problem with Cayley parameterization,” BMVC 2015.
【非特許文献４】Y. Zheng, et al., “A General and Simple Method for Camera Pose and Focal Length Determination,” CVPR 2014.
【非特許文献５】M. Bujnak, et al., “New efficient solution to the absolute pose problem for camera with unknown focal length and radial distortion,” ACCV 2010.
【非特許文献６】Z. Kukelova, et al., “Real-time solution to the absolute pose problem with unknown radial distortion and focal length,” ICCV 2013.
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、上述の非特許文献１〜６に開示された方法には、以下のような問題がある。まず、非特許文献１および２に開示された方法では、安定的な計算には数十枚の画像が必要であり、また、非線形最適化の際に望ましくない局所解に収束することがあり、内部パラメータおよび外部パラメータの推定精度が低いという問題がある。
【０００９】
また、非特許文献３に開示された方法では、固定焦点距離のカメラが想定されており、ズームレンズを搭載したカメラが用いられた場合に、内部パラメータおよび外部パラメータを推定できないという問題がある。一方、非特許文献４に開示された方法では、ズームレンズを搭載したカメラのズームを変化させた際の焦点距離は推定できるものの、ズームの変化と同時に生じるレンズ歪み係数の変化に対応できない。よって、非特許文献４に開示された方法にも、内部パラメータおよび外部パラメータの推定精度が低いという問題がある。
【００１０】
更に、非特許文献５または非特許文献６に開示された方法では、一般に複数の解が得られ、一意に解を選択できず、加えて、４組以上または５組以上の対応点を入力として最小二乗的な意味で最適な解を推定することもできない。従って、非特許文献５または非特許文献６に開示された方法にも、内部パラメータおよび外部パラメータの推定精度が低いという問題がある。
【００１１】
このように、上述した非特許文献１〜５に開示された方法には、入力である複数の対応点から、レンズ歪み係数を未知として含む内部パラメータと外部パラメータとを高精度に推定することは難しいという問題がある。
【００１２】
本発明の目的の一例は、上記問題を解消し、レンズ歪み係数を含む内部パラメータと外部パラメータとを高精度に推定し得る、カメラパラメータ推定装置、カメラパラメータ推定方法、およびプログラムを提供することにある。
【課題を解決するための手段】
【００１３】
上記目的を達成するため、本発明の一側面におけるカメラパラメータ推定装置は、カメラによる対象物の撮影時におけるカメラパラメータを推定するための装置であって、
前記対象物に関する５組以上の３次元点と、前記３次元点それぞれに対応する画像上の２次元点との変換を表す誤差関数を、前記カメラパラメータの依存関係に基づいて、複数の部分問題に分割し、前記複数の部分問題それぞれにおいて停留点を算出する、パラメータ算出部と、
前記誤差関数によって得られる誤差を最小とする停留点を求め、求めた前記停留点を最適なカメラパラメータとして出力する、パラメータ最適化部と、
を備えている、ことを特徴とする。
【００１４】
また、上記目的を達成するため、本発明の一側面におけるカメラパラメータ推定方法は、カメラによる対象物の撮影時におけるカメラパラメータを推定するための方法であって、
（ａ）前記対象物に関する５組以上の３次元点と、前記３次元点それぞれに対応する画像上の２次元点との変換を表す誤差関数を、前記カメラパラメータの依存関係に基づいて、複数の部分問題に分割し、前記複数の部分問題それぞれにおいて停留点を算出する、ステップと、
（ｂ）前記誤差関数によって得られる誤差を最小とする停留点を求め、求めた前記停留点を最適なカメラパラメータとして出力する、ステップと、
を有する、ことを特徴とする。
【００１５】
更に、上記目的を達成するため、本発明の一側面におけるプログラムは、コンピュータによって、カメラによる対象物の撮影時におけるカメラパラメータを推定するためのプログラムであって、
前記コンピュータに、
（ａ）前記対象物に関する５組以上の３次元点と、前記３次元点それぞれに対応する画像上の２次元点との変換を表す誤差関数を、前記カメラパラメータの依存関係に基づいて、複数の部分問題に分割し、前記複数の部分問題それぞれにおいて停留点を算出する、ステップと、
（ｂ）前記誤差関数によって得られる誤差を最小とする停留点を求め、求めた前記停留点を最適なカメラパラメータとして出力する、ステップと、
を実行させることを特徴とする。
【発明の効果】
【００１６】
以上のように、本発明によれば、レンズ歪み係数を含む内部パラメータと外部パラメータとを高精度に推定することができる。
【図面の簡単な説明】
【００１７】
【図１】図１は、本発明の実施の形態におけるカメラパラメータ推定装置の概略構成を示すブロック図である。
【図２】図２は、本発明の実施の形態におけるカメラパラメータ推定装置の具体的構成を示すブロック図である。
【図３】図３は、本発明の実施の形態におけるカメラパラメータ推定装置１０の動作を示すフロー図である。
【図４】図４は、本発明の実施の形態におけるカメラパラメータ推定装置１０を実現するコンピュータの一例を示すブロック図である。
【発明を実施するための形態】
【００１８】
（実施の形態）
以下、本発明の実施の形態における、カメラパラメータ推定装置、カメラパラメータ推定方法、及びプログラムについて、図１〜図４を参照しながら説明する。
【００１９】
［装置構成］
最初に、本発明の実施の形態におけるカメラパラメータ推定装置の概略構成及び機能について図１を用いて説明する。図１は、本発明の実施の形態におけるカメラパラメータ推定装置の概略構成を示すブロック図である。
【００２０】
図１に示す本実施の形態におけるカメラパラメータ推定装置１０は、カメラによる対象物の撮影時におけるカメラパラメータを推定する装置である。図１に示すように、カメラパラメータ推定装置１０は、パラメータ算出部２０と、パラメータ最適化部３０とを備えている。
【００２１】
パラメータ算出部２０は、まず、対象物に関する５組以上の３次元点と、３次元点それぞれに対応する画像上の２次元点との変換を表す誤差関数を、カメラパラメータの依存関係に基づいて、複数の部分問題に分割する。次いで、パラメータ算出部２０は、複数の部分問題それぞれにおいて停留点を算出する。
【００２２】
パラメータ最適化部３０は、まず、誤差関数によって得られる誤差を最小とする停留点を求める。次いで、パラメータ最適化部３０は、求めた停留点を最適なカメラパラメータとして出力する。
【００２３】
このように、本実施の形態では、誤差関数は複数の部分問題に分割され、部分問題毎に最適解が計算される。よって、本実施の形態によれば、カメラパラメータ、即ち、レンズ歪み係数を含む内部パラメータと外部パラメータとが高精度に推定される。
【００２４】
続いて、図２を用いて、本発明の実施の形態におけるカメラパラメータ推定装置の構成及び機能について、更に具体的に説明する。図２は、本発明の実施の形態におけるカメラパラメータ推定装置の具体的構成を示すブロック図である。
【００２５】
図２に示すように、本実施の形態では、カメラパラメータ推定装置１０には、撮影の対象物に関する５組以上の３次元点と、３次元点それぞれに対応する画像上の２次元点とが、上述した対応点として入力される。
【００２６】
本実施の形態において、カメラパラメータには、カメラの外部に関する外部パラメータと、カメラの内部に関する内部パラメータとが含まれる。具体的には、外部パラメータは、従来と同様に、カメラの３次元区間における位置及び姿勢のいずれか一方または両方を含む。内部パラメータは、カメラにおける、レンズの焦点距離、光学中心、アスペクト比、せん断係数、及びレンズ歪み係数のうち、少なくとも１つを含む。なお、本実施の形態では、カメラパラメータには、外部パラメータのみが含まれていても良い。
【００２７】
本実施の形態では、カメラパラメータ推定装置１０は、５組以上の対応点に基づいて、カメラの外部パラメータ、更にはカメラの外部パラメータに加えて内部パラメータも推定する情報処理装置である。また、本実施の形態では、図２に示すように、パラメータ算出部２０には、分割処理部２１と、部分パラメータ算出部２２と、残存パラメータ算出部２３とが備えられている。
【００２８】
分割処理部２１は、対応点に含まれる３次元点と２次元点との所定の変換を表す誤差関数を、カメラパラメータの依存関係に基づいて、複数の部分問題に分割する。ここで「カメラパラメータ依存関係に基づく分割」とは、誤差関数を構成する複数の項のうち、例えば外部パラメータのみが変数として現れる項、または、内部パラメータのみが変数として現れる項のように、未知のパラメータ同士の影響を受けない、もしくは影響が少ない項を選出することを表す。
【００２９】
部分パラメータ算出部２２は、対応点を入力として、複数の部分問題のうち、外部パラメータの一部のみに依存する１つの部分問題について停留点を算出する。
【００３０】
残存パラメータ算出部２３は、部分パラメータ算出部２２で算出された停留点と、対応点とを入力として、これらに基づいて、停留点が未だ算出されていない部分問題について停留点を算出する。
【００３１】
パラメータ最適化部３０は、本実施の形態では、部分パラメータ算出部２２と残存パラメータ算出部２３とで算出された停留点と対応点とを入力として、これらに基づいて、分割される前の元の誤差関数を最適化することによって、誤差を最小とする停留点を求める。そして、パラメータ最適化部３０は、求めた停留点を最適なカメラパラメータとして出力する。
【００３２】
［装置動作］
次に、本発明の実施の形態におけるカメラパラメータ推定装置１０の動作について図３を用いて説明する。図３は、本発明の実施の形態におけるカメラパラメータ推定装置１０の動作を示すフロー図である。以下の説明においては、適宜図１を参酌する。また、本実施の形態では、カメラパラメータ推定装置１０を動作させることによって、カメラパラメータ推定方法が実施される。よって、本実施の形態におけるカメラパラメータ推定方法の説明は、以下のパラメータ推定装置１０の動作説明に代える。
【００３３】
図３に示すように、まず、カメラパラメータ推定装置１０に５組以上の対応点が入力されると、分割処理部２１は、対応点に含まれる３次元点と２次元点との所定の変換を表す誤差関数を、カメラパラメータの依存関係に基づいて、複数の部分問題に分割する（ステップＳ１１）。
【００３４】
本実施の形態において、入力される対応点は、推定すべきカメラパラメータにおける、いわゆる退化配置（Degenerate ConfigurationまたはCritical Configurationとも呼ばれる）ではないとする。なぜならば、そのような対応点に対しては理論的にカメラパラメータを推定できないためである。また、座標値に関して、３次元点も２次元点も共にすべて異なるとする。これは、例えば、１つの３次元点が複数の異なる２次元点に対応することは現実的に起こりえないためである。カメラパラメータ推定装置１０において、これらの不適切な対応点は以降の処理から除外される。また、本実施の形態は、カメラパラメータ推定装置１０の代わりに、対応点を外部からカメラパラメータ推定装置に入力する別の装置（図示せず）が、対応点のエラー判定等を行う態様であっても良い。
【００３５】
次に、部分パラメータ算出部２２は、ステップＳ１１の分割によって得られた部分問題の中で外部パラメータの一部のみに依存する１つの部分問題について、すべての停留点を算出する（ステップＳ１２）。
【００３６】
ここで、停留点とは、最適化問題における目的関数（より正確にはラグランジェ関数であるが、以下単に「目的関数」と表記する）の勾配ベクトルがゼロとなる未知数の点のことを意味する。すなわち、停留点は、勾配ベクトルを連立方程式とみなすと、その連立方程式の解のことである。ステップＳ１２においては、目的関数は外部パラメータの一部のみに依存する１つの部分問題を指し、停留点は当該部分問題の勾配ベクトルをゼロとする外部パラメータを指す。
【００３７】
一般的に、連立方程式の解は複素数になることもありえるが、すべてのカメラパラメータは必ず実数であるため、本実施の形態においては、実数解のみが対象となる。また、カメラパラメータ推定装置１０は、各ステップの計算結果が複素数のみとなった場合には、解なしのフラグを出して以降の処理を終了してもよい。例えば、誤った３次元点と２次元点との組が対応点に含まれていたり、数値計算の精度が不足して実数解が複素数解として計算されたりする場合に、連立方程式の解が複素数になることが起こりうる。
【００３８】
次いで、残存パラメータ算出部２３は、対応点と、ステップＳ１２で算出された停留点とを入力として、残りのカメラパラメータを未知数とするその他の部分問題、即ち、停留点が未だ算出されていない部分問題における、すべての停留点を算出する（ステップＳ１３）。
【００３９】
ここで、残りのカメラパラメータとは、推定したい全カメラパラメータから、ステップＳ１２で停留点が算出された部分問題に依存する、外部パラメータの一部を除いた、カメラパラメータのことである。元の分割される前の誤差関数においては、すべてのカメラパラメータが未知数であるが、ステップＳ１３における部分問題においては、ステップＳ１２で算出された停留点（カメラパラメータ）が既知の定数として扱われる。すなわち、ステップＳ１３では、ステップＳ１２で算出されなかった残りのカメラパラメータが、ステップＳ１２では対象とならなかった部分問題の停留点として算出される。
【００４０】
また、ステップＳ１２で出力され停留点が複数の場合、残存パラメータ算出部２３は、複数の停留点それぞれ毎に、その他の部分問題の停留点の算出を行う。すなわち、例えば、ステップＳ１２でＫ個の停留点が算出された場合は、ステップＳ１３は、その他の部分問題をＫ回解くことになる。
【００４１】
最後に、パラメータ最適化部３０は、ステップＳ１１とステップＳ１２とで出力されたすべての停留点と対応点とを入力として、元の誤差関数の最適化を行い、誤差が最小となる解（停留点）を求め、求めた解をカメラパラメータとして出力する（ステップＳ１３）。
【００４２】
ステップＳ１３では、パラメータ最適化部３０は、更に、入力または出力されたカメラパラメータが物理的または理論的に意味のあるかどうかの判定を行うこともできる。物理的な意味とは、例えば、すべての３次元点はカメラの前面に存在しなくてはならないことをいう。従って、パラメータ最適化部３０は、一部またはすべての３次元点がカメラの裏側に存在するような不適切なカメラパラメータを除外してもよい。また、一部のカメラパラメータが取りうる値の範囲が予めわかっている場合、パラメータ最適化部３０は、それを閾値として不適切な入出力を除外してもよい。例えば、ズームレンズの仕様がわかっていれば、推定される焦点距離の上限と下限を事前に決めることができるので、パラメータ最適化部３０は、この上限及び下限から外れる入出力を除外する。
【００４３】
理論的な意味とは、例えば、計算結果が誤差関数の２次最適性を満たさないことをいう。パラメータ最適化部３０は、計算されたカメラパラメータが誤差関数の２次最適性を満たさない場合は、このカメラパラメータを除外する。
【００４４】
このように不適切なカメラパラメータが除外された場合、パラメータ最適化部３０は、物理的または理論的に意味のある値の範囲内での、誤差最小となる解を出力する。
【００４５】
もし上記の条件を満たす解が複数存在し、それらの誤差が同一とみなせる場合には、パラメータ最適化部３０は、すべての解を出力してもよい。この場合、大域的最適解が複数存在するとみなせるため、出力されるカメラパラメータも複数である。
【００４６】
［具体例］
ここで、本実施の形態の具体例について以下に説明する。なお、以下の具体例では、外部パラメータは、それぞれ３自由度の位置と姿勢とする。また、ズームレンズにおいてズームを変更した際に影響を受けるのは焦点距離とレンズ歪み係数であるため、内部パラメータは１自由度の焦点距離と３自由度のレンズ歪み係数とする。
【００４７】
レンズ歪みの種類は樽型歪（radial distortion）とし、計算モデルはAndrew Fitzgibbonが提案したdivision modelとする。それ以外の内部パラメータは、光学中心については画像の中心（つまり画像の大きさの半分）と仮定し、アスペクト比については１と仮定し、せん断係数については０と仮定する。これらの仮定は、上述した非特許文献４、５、６でも採用されている。昨今のデジタルカメラは、高精度に製造されているため、上記の仮定は妥当な仮定である。
【００４８】
但し、どのカメラパラメータを未知または既知とするかはこの限りではなく、利用シーンに応じて当業者が自由に変更することができる。また、レンズ歪みの種類および計算モデルについても同様である。
【００４９】
また、以下の説明では、上添字の「Ｔ」は行列およびベクトルの転置を表し、「０」はゼロ行列およびゼロベクトルを表し、「Ｉ」は単位行列を表し、「||  ||」はベクトルのＬ２ノルムを表す。外部パラメータにおいては、カメラ姿勢を３×３の回転行列Ｒで表し、カメラ位置を３×１のベクトルｔで表す。内部パラメータにおいては、焦点距離を「ｆ」で表し、焦点距離の逆数１／ｆを「ｗ」で表記し、レンズ歪み係数を３×１のベクトル「ｋ」で表す。また、ベクトルｋを下記の数１で表す。
【００５０】
【数１】
2017163386.tif 000003
【００５１】
対応点は、入力点数をｎ組（ｎ≧５）とし、i番目の３次元座標を「」、division modelにより斉次化した２次元座標を「ｍｉ」と表記する。２次元座標ｍｉは下記の数２で表される。
【００５２】
【数２】
2017163386.tif 000004
【００５３】
また、誤差関数の定義を行う。誤差関数は、Ｘｉとのｍｉとの変換関係を表しており、以下の数３で定義される。
【００５４】
【数３】
2017163386.tif 000005
【００５５】
上記数３において、「〜」は定数倍に不定なことを表している。「Ｋ」は、［１，１，ｗ］を対角成分とする３×３の行列である。
【００５６】
また、上記数３より、ｎ組の対応点に対する誤差関数は以下の数４で定義される。
【００５７】
【数４】
2017163386.tif 000006
【００５８】
ここで、［ｍｉ］ｘはベクトルの外積を示す３×３の歪対称行列であり、具体的な成分は下記の数５で表される。
【００５９】
【数５】
2017163386.tif 000007
【００６０】
上記数４は、代数学的誤差と呼ばれ、上記数３の両辺のベクトル方向が互いに近しくなることを要求している。
【００６１】
次に、部分問題の定義を行う。［ｍｉ］ｘの行ベクトルを下記の数６に示す通りとすると上記数４は以下の数９と等価である。なお、制約条件の記述は省略している。
【００６２】
【数６】
2017163386.tif 000008
【００６３】
【数７】
2017163386.tif 000009
【００６４】
【数８】
2017163386.tif 000010
【００６５】
【数９】
2017163386.tif 000011
【００６６】
ここで、上記数９の各項を独立の最適化問題とみなすと、以下の数１０〜数１２（制約条件の記述は省略）で表される３つの部分問題が得られる。
【００６７】
【数１０】
2017163386.tif 000012
【００６８】
【数１１】
2017163386.tif 000013
【００６９】
【数１２】
2017163386.tif 000014
【００７０】
次に、各部分問題の具体的な解法について説明する。まず、上記数１２に着目すると、下記の数１３が成立するため、上記数１２を下記の数１４で表すことができる。
【００７１】
【数１３】
2017163386.tif 000015
【００７２】
【数１４】
2017163386.tif 000016
【００７３】
上記数１４において、ｒｊＴはＲの第j行、ｔｊはｔの第j成分である。また、上記数１４は、内部パラメータを含まず、ｔ３を除く外部パラメータのみで表現されている。また、Ｒに関する制約条件を考慮すると、｛ｒ１，ｒ２，ｔ１，ｔ２｝の自由度は５であるため、対応点数ｎが５組以上であれば、上記数１２を解くことができる。
【００７４】
また、上記数１４は、内部パラメータを既知とする非特許文献３に開示されている形式と同様の形式であるため、非特許文献３で用いられている様々な多項式表現による方法により、上記数１２を解くことができる。すなわち、Ｒを４変数の正規化四元数、３変数の非正規化四元数（Cayley parameterization）、９変数の回転行列などで表現し、４０×４０の行列固有値問題へ帰着させる方法によれば、上記数１２において、最大で４０個の実数解を得ることができる。
【００７５】
また、上記数３において｛ｒ１，ｒ２，ｗ｝も｛−ｒ１，−ｒ２，−ｗ｝も解になる。このことから、｛ｒ１，ｒ２｝に関する正負符号を固定する制約条件（例えば、ｒ１のうち１つを正の定数とおく）を導入すれば、符号が負の場合でも正の場合と同様の方法で２０×２０の行列固有値問題へ帰着させることもできる。この場合、ｒ３はｒ１とｒ２との外積として、ｒ１とｒ２との推定後に算出すればよい。連立多項式問題から行列固有値問題を導く方法は非特許文献３に詳しく記載されているため、ここでは記述を省略する。
【００７６】
上記数１２が解けると、ｔ３を除く外部パラメータが求まっているため、残るカメラパラメータはｔ３，ｗ，ｋである。そこで、ｔ３を除く外部パラメータを既知とすると、上記数１０は下記の数１５で表され、数１１は下記の数１６で表される。
【００７７】
【数１５】
2017163386.tif 000017
【００７８】
【数１６】
2017163386.tif 000018
【００７９】
対応点数ｎは５以上であり、また、残るパラメータ数はｔ３，ｗ，ｋ＝［ｋ１，ｋ２，ｋ３］の３つであるから、上記数１５及び数１６より、ｔ３，ｗ，ｋを計算できる。
【００８０】
例として、上記数１５を用いる場合について説明する。上記数１５を展開すると、以下の数１７が得られる。
【００８１】
【数１７】
2017163386.tif 000019
【００８２】
上記数１７は、［ｔ３，ｗ，ｋ］Ｔに関する線形最小二乗問題とみなされるため、上記数１７を解くことは容易である。
【００８３】
このように、すべてのカメラパラメータが推定できたが、これらは上記数１０〜数１２で表される部分問題の解である。そこで、これらを初期値としてニュートン法などの非線形最適化法を実行して、上記数４または数９を解けば、最終的に最適なカメラパラメータが推定できる。
【００８４】
ここで、上述の具体例を図２に示したステップＳ１１〜Ｓ１４に沿って以下に説明する。
【００８５】
［ステップＳ１１］
まず、ステップＳ１１において、分割処理部２１は、上記数９で表される誤差関数を分割し、上記数１０、数１１、及び数１２の部分問題とする。
【００８６】
［ステップＳ１２］
次に、ステップＳ１２において、部分パラメータ算出部２２は、上記数１２で表される部分問題を、連立多項式問題として解いて、その停留点を出力する。連立多項式の解法としては、非特許文献３に記されている方法に加えて、Ｆ４アルゴリズムなどが利用されてもよい。また、このとき、部分パラメータ算出部２２は、すべての停留点を出力してもよいし、上記数１２を最小化する唯一の停留点を出力してもよい。
【００８７】
［ステップＳ１３］
次に、ステップＳ１３において、残存パラメータ算出部２３は、部分パラメータ算出部２２により計算された外部パラメータの一部を表す停留点を入力として、上記数１５および上記数１６を線形最小二乗法により解いて、得られた残りの停留点、すなわち、残りの外部パラメータと内部パラメータとを出力する。ステップＳ１２で算出された停留点が複数であり、複数の停留点が入力された場合は、残存パラメータ算出部２３は、入力された停留点毎に、上記数１５及び上記数１６を解く。このとき、残存パラメータ算出部２３は、すべてのパラメータを出力してもよいし、上記数４、数１５、及び数１６のいずれかを最小化する唯一の停留点を出力してもよい。
【００８８】
［ステップＳ１４］
最後に、ステップＳ１４において、パラメータ最適化部３０は、部分パラメータ算出部２２と残存パラメータ算出部２３とで出力された１組以上のカメラパラメータを初期値として、上記数４の非線形最適化により誤差が最小となる最適解を出力する。ステップＳ１４の実行後、カメラパラメータ推定装置１０は動作を停止する。
【００８９】
また、ステップＳ１４において、パラメータ最適化部３０は、最適化前に入力パラメータを用いて、上記数４の誤差を計算し、誤差が閾値以下であれば最適化をせずに、入力パラメータをそのまま出力してもよい。また、ステップＳ１４においては、誤差関数として、上記数４以外にも、物理的な距離を表すいわゆる再投影誤差関数が採用されていてもよい。この場合、上記数４を最適化する場合よりも計算時間はかかるが、より高精度な解が得られる可能性がある。
【００９０】
［実施の形態における効果］
以上のとおり、本実施形態によれば、レンズ歪み係数を含む内部パラメータと外部パラメータとを高精度に推定可能である。その理由は、以下のとおりである。
【００９１】
上記数４および数９は２つのベクトルｍｉとＫ（ＲＸｉ＋ｔ）との方向が近しくなることを表している。換言すれば、上記数４および数９の最適解が求まったとすると、［ｍｉ］ｘＫ（ＲＸｉ＋ｔ）で表されるベクトルの各成分はそれぞれ０に近くなると期待される。そこで、上記数９を直接に解く代わりに、上記数９の各項を独立の制約付き最適化問題（上記数１０〜数１２の部分問題）とみなして解いたとしても、それぞれの最適解は原問題である上記数９の最適解に非常に近いと考えられる。更に、通常であれば計算量は未知数の数に依存して増加するが、本実施の形態では、外部パラメータのみに依存する部分問題（上記数１２）から順に解くことで、従来方法と同等以下の低演算量で実行が可能である。以上の理由により、本実施の形態によれば、内部パラメータと外部パラメータとが高精度に推定される。
【００９２】
［変形例］
本実施の形態は、上述した例に限定されるものではない。本実施の形態は、上述した例に対して、いわゆる当業者が理解し得る多様な変更を適用することが可能である。例えば、本実施の形態は、以下の変形例に示す形態によっても実施可能である。
【００９３】
変形例１
本実施の形態において、未知とする内部パラメータは、特に限定されるものではない。例えば、非特許文献４に開示されているように、推定すべき内部パラメータが焦点距離のみの場合、上記数１７においてｋ＝０として解けばよい。また、非特許文献３に開示されているように、すべての内部パラメータが既知である場合、上記数１７において、ｋ＝０、ｗ＝１／ｆとして解けばよい。変形例１においても、未知となる内部パラメータをどのように設定しようとも、上記数１２を用いて外部パラメータの一部を算出してから、上記数１７を解いて、その他のパラメータを推定するという方法は変更する必要はない。
【００９４】
変形例２
本実施の形態では、カメラパラメータ推定装置１０は、一部の外部パラメータを外部より受け付けるために、外部パラメータ入力部を備えていてもよい。外部パラメータ入力部は、例えば、外部の加速度センサから、重力方向における加速度を特定するセンサ値を受け付けたり、ジャイロセンサから、角速度を特定するセンサ値を受け付けたりしてもよい。後者の場合、外部パラメータ入力部は、回転行列を算出し、得られた回転行列をパラメータ算出部２０に入力する。変形例２においては、推定すべき外部パラメータ数は減少する。但し、変形例２においても、上記数１２で未知の外部パラメータを算出してから、上記数１７を解いて、内部パラメータを推定するという方法は変更する必要はない。
【００９５】
変形例３
本実施の形態において、外部パラメータはカメラの位置姿勢に限定されない。例えば、静止したカメラの前で動く対象物を撮影した場合、外部パラメータとして、対象物の位置姿勢を表すパラメータが推定される。カメラが動くか、対象物が動くかは相対的に等価なため、解くべき問題は変わらない。
【００９６】
変形例４
本実施の形態において、撮影に用いられるカメラは単眼カメラに限定されない。上記数３では、Ｘｉとのｍｉとの変換関係を、いわゆる単眼のピンホールカメラモデルを用いて説明しているが、本実施の形態では、複数の焦点を持つ一般カメラモデルが用いられていてもよい。焦点位置が事前に既知であれば、上述した例と同様の方法でカメラパラメータを推定可能である。
【００９７】
［プログラム］
本実施の形態におけるプログラムは、コンピュータに、図３に示すステップＳ１１〜Ｓ１４を実行させるプログラムであれば良い。このプログラムをコンピュータにインストールし、実行することによって、本実施の形態におけるカメラパラメータ推定装置１０とカメラパラメータ推定方法とを実現することができる。この場合、コンピュータのＣＰＵ（Central Processing Unit）は、パラメータ算出部２０及びパラメータ最適化部３０部として機能し、処理を行なう。
【００９８】
また、本実施の形態におけるプログラムは、複数のコンピュータによって構築されたコンピュータシステムによって実行されても良い。この場合は、例えば、各コンピュータが、それぞれ、パラメータ算出部２０、及びパラメータ最適化部３０のいずれかとして機能しても良い。
【００９９】
［物理構成］
ここで、本実施の形態におけるプログラムを実行することによって、カメラパラメータ推定装置１０を実現するコンピュータについて図４を用いて説明する。図４は、本発明の実施の形態におけるカメラパラメータ推定装置１０を実現するコンピュータの一例を示すブロック図である。
【０１００】
図４に示すように、コンピュータ１１０は、ＣＰＵ１１１と、メインメモリ１１２と、記憶装置１１３と、入力インターフェイス１１４と、表示コントローラ１１５と、データリーダ／ライタ１１６と、通信インターフェイス１１７とを備える。これらの各部は、バス１２１を介して、互いにデータ通信可能に接続される。
【０１０１】
ＣＰＵ１１１は、記憶装置１１３に格納された、本実施の形態におけるプログラム（コード）をメインメモリ１１２に展開し、これらを所定順序で実行することにより、各種の演算を実施する。メインメモリ１１２は、典型的には、ＤＲＡＭ（Dynamic Random Access Memory）等の揮発性の記憶装置である。また、本実施の形態におけるプログラムは、コンピュータ読み取り可能な記録媒体１２０に格納された状態で提供される。なお、本実施の形態におけるプログラムは、通信インターフェイス１１７を介して接続されたインターネット上で流通するものであっても良い。
【０１０２】
また、記憶装置１１３の具体例としては、ハードディスクドライブの他、フラッシュメモリ等の半導体記憶装置が挙げられる。入力インターフェイス１１４は、ＣＰＵ１１１と、キーボード及びマウスといった入力機器１１８との間のデータ伝送を仲介する。表示コントローラ１１５は、ディスプレイ装置１１９と接続され、ディスプレイ装置１１９での表示を制御する。
【０１０３】
データリーダ／ライタ１１６は、ＣＰＵ１１１と記録媒体１２０との間のデータ伝送を仲介し、記録媒体１２０からのプログラムの読み出し、及びコンピュータ１１０における処理結果の記録媒体１２０への書き込みを実行する。通信インターフェイス１１７は、ＣＰＵ１１１と、他のコンピュータとの間のデータ伝送を仲介する。
【０１０４】
また、記録媒体１２０の具体例としては、ＣＦ（Compact Flash（登録商標））及びＳＤ（Secure Digital）等の汎用的な半導体記憶デバイス、フレキシブルディスク（Flexible Disk）等の磁気記憶媒体、又はＣＤ−ＲＯＭ（Compact Disk Read Only Memory）などの光学記憶媒体が挙げられる。
【０１０５】
なお、本実施の形態におけるカメラパラメータ推定装置１０は、プログラムがインストールされたコンピュータではなく、各部に対応したハードウェアを用いることによっても実現可能である。更に、カメラパラメータ推定装置１０は、一部がプログラムで実現され、残りの部分がハードウェアで実現されていてもよい。また、カメラパラメータ推定装置１０は、各部が、それぞれ別々のユニットとして実現されていてもよい。
【０１０６】
上述した実施の形態の一部又は全部は、以下に記載する（付記１）〜（付記１２）によって表現することができるが、以下の記載に限定されるものではない。
【０１０７】
（付記１）
カメラによる対象物の撮影時におけるカメラパラメータを推定するための装置であって、
前記対象物に関する５組以上の３次元点と、前記３次元点それぞれに対応する画像上の２次元点との変換を表す誤差関数を、前記カメラパラメータの依存関係に基づいて、複数の部分問題に分割し、前記複数の部分問題それぞれにおいて停留点を算出する、パラメータ算出部と、
前記誤差関数によって得られる誤差を最小とする停留点を求め、求めた前記停留点を最適なカメラパラメータとして出力する、パラメータ最適化部と、
を備えている、ことを特徴とするカメラパラメータ推定装置。
【０１０８】
（付記２）
前記カメラパラメータが、前記カメラの外部に関する外部パラメータを含み、
前記パラメータ算出部が、
前記複数の部分問題のうち、前記外部パラメータの一部のみに依存する１つの部分問題について、前記停留点を算出し、
更に、算出した前記停留点と、前記３次元点及び前記２次元点とに基づいて、前記停留点が未だ算出されていない前記部分問題について、前記停留点を算出する、
付記１に記載のカメラパラメータ推定装置。
【０１０９】
（付記３）
前記パラメータ最適化部が、前記パラメータ算出部によって算出された前記停留点と、前記３次元点及び前記２次元点とに基づいて、前記誤差関数を最適化することによって、前記誤差を最小とする停留点を求める、
付記１または２に記載のカメラパラメータ推定装置。
【０１１０】
（付記４）
前記外部パラメータが、前記カメラの３次元区間における位置及び姿勢のいずれか一方または両方を含み、
前記カメラパラメータが、前記外部パラメータに加えて、前記カメラの内部に関する内部パラメータを含み、前記内部パラメータは、前記カメラにおける、レンズの焦点距離、光学中心、アスペクト比、せん断係数、及びレンズ歪み係数のうち、少なくとも１つを含む、
付記２に記載のカメラパラメータ推定装置。
【０１１１】
（付記５）
カメラによる対象物の撮影時におけるカメラパラメータを推定するための方法であって、
（ａ）前記対象物に関する５組以上の３次元点と、前記３次元点それぞれに対応する画像上の２次元点との変換を表す誤差関数を、前記カメラパラメータの依存関係に基づいて、複数の部分問題に分割し、前記複数の部分問題それぞれにおいて停留点を算出する、ステップと、
（ｂ）前記誤差関数によって得られる誤差を最小とする停留点を求め、求めた前記停留点を最適なカメラパラメータとして出力する、ステップと、
を有する、ことを特徴とするカメラパラメータ推定方法。
【０１１２】
（付記６）
前記カメラパラメータが、前記カメラの外部に関する外部パラメータを含み、
前記（ａ）のステップにおいて、
前記複数の部分問題のうち、前記外部パラメータの一部のみに依存する１つの部分問題について、前記停留点を算出し、
更に、算出した前記停留点と、前記３次元点及び前記２次元点とに基づいて、前記停留点が未だ算出されていない前記部分問題について、前記停留点を算出する、
付記５に記載のカメラパラメータ推定方法。
【０１１３】
（付記７）
前記（ｂ）のステップにおいて、前記パラメータ算出部によって算出された前記停留点と、前記３次元点及び前記２次元点とに基づいて、前記誤差関数を最適化することによって、前記誤差を最小とする停留点を求める、
付記５または６に記載のカメラパラメータ推定方法。
【０１１４】
（付記８）
前記外部パラメータが、前記カメラの３次元区間における位置及び姿勢のいずれか一方または両方を含み、
前記カメラパラメータが、前記外部パラメータに加えて、前記カメラの内部に関する内部パラメータを含み、前記内部パラメータは、前記カメラにおける、レンズの焦点距離、光学中心、アスペクト比、せん断係数、及びレンズ歪み係数のうち、少なくとも１つを含む、
付記６に記載のカメラパラメータ推定方法。
【０１１５】
（付記９）
コンピュータによって、カメラによる対象物の撮影時におけるカメラパラメータを推定するためのプログラムであって、
前記コンピュータに、
（ａ）前記対象物に関する５組以上の３次元点と、前記３次元点それぞれに対応する画像上の２次元点との変換を表す誤差関数を、前記カメラパラメータの依存関係に基づいて、複数の部分問題に分割し、前記複数の部分問題それぞれにおいて停留点を算出する、ステップと、
（ｂ）前記誤差関数によって得られる誤差を最小とする停留点を求め、求めた前記停留点を最適なカメラパラメータとして出力する、ステップと、
を実行させるプログラム。
【０１１６】
（付記１０）
前記カメラパラメータが、前記カメラの外部に関する外部パラメータを含み、
前記（ａ）のステップにおいて、
前記複数の部分問題のうち、前記外部パラメータの一部のみに依存する１つの部分問題について、前記停留点を算出し、
更に、算出した前記停留点と、前記３次元点及び前記２次元点とに基づいて、前記停留点が未だ算出されていない前記部分問題について、前記停留点を算出する、
付記９に記載のプログラム。
【０１１７】
（付記１１）
前記（ｂ）のステップにおいて、前記パラメータ算出部によって算出された前記停留点と、前記３次元点及び前記２次元点とに基づいて、前記誤差関数を最適化することによって、前記誤差を最小とする停留点を求める、
付記９または１０に記載のプログラム。
【０１１８】
（付記１２）
前記外部パラメータが、前記カメラの３次元区間における位置及び姿勢のいずれか一方または両方を含み、
前記カメラパラメータが、前記外部パラメータに加えて、前記カメラの内部に関する内部パラメータを含み、前記内部パラメータは、前記カメラにおける、レンズの焦点距離、光学中心、アスペクト比、せん断係数、及びレンズ歪み係数のうち、少なくとも１つを含む、
付記１０に記載のプログラム。
【産業上の利用可能性】
【０１１９】
以上のように、本発明によれば、レンズ歪み係数を含む内部パラメータと外部パラメータとを高精度に推定することができる。本発明は、カメラパラメータの推定が求められる分野、例えば、ロボットの制御、画像処理等に有効である。
【符号の説明】
【０１２０】
１０  カメラパラメータ推定装置
２０  パラメータ算出部
２１  分割処理部
２２  部分パラメータ算出部
２２  残存パラメータ算出部
３０  パラメータ最適化部
１１０  コンピュータ
１１１  ＣＰＵ
１１２  メインメモリ
１１３  記憶装置
１１４  入力インターフェイス
１１５  表示コントローラ
１１６  データリーダ／ライタ
１１７  通信インターフェイス
１１８  入力機器
１１９  ディスプレイ装置
１２０  記録媒体
１２１  バス
【図１】
2017163386.tif 000020
【図２】
2017163386.tif 000021
【図３】
2017163386.tif 000022
【図４】
2017163386.tif 000023