JP6525195B2 - 魚眼カメラの内部パラメータの、推定方法、推定装置、及び推定プログラム - Google Patents

Description

本発明は、魚眼レンズを備えた魚眼カメラの内部パラメータの、推定方法、推定装置、及び推定プログラムに関するものである。

従来、魚眼カメラを用いて、例えばストライプを含むパターン等を有する被写体を撮像し、撮像により得られた画像に基づいて、魚眼カメラの内部パラメータを推定する技術が、知られている（例えば、特許文献１）。

特開2007-192832号公報

しかしながら、上述した技術では、上記のように特別な構成を持った被写体を事前に用意する必要があり、また、そのような被写体を魚眼カメラで撮像するにあたって、魚眼カメラに対する被写体の向きや距離等に応じて、内部パラメータの推定結果がばらつくおそれがあった。

本発明は、上述した課題を解決するためのものであり、簡単かつ高精度に、魚眼カメラの内部パラメータを推定できる、魚眼カメラの内部パラメータの、推定方法、推定装置、及び推定プログラムを提供することを目的とするものである。

本発明の推定方法は、魚眼レンズを備えた魚眼カメラの内部パラメータの推定方法であって、
前記魚眼カメラが、複数の視覚的な特徴点を含む空間内で、該魚眼カメラの内部を通るとともに該魚眼カメラの光軸に対してほぼ垂直に延在する回転軸線の周りで回転される間に、前記魚眼カメラによって、複数の異なるタイミングのそれぞれで前記複数の特徴点を撮像することにより、複数の魚眼画像を得る、撮像ステップと、
前記各特徴点に対応する前記各魚眼画像中の観測点の集合から得られる軌跡を含む、軌跡画像を生成する、軌跡画像生成ステップと、
前記軌跡画像に基づいて、前記魚眼カメラの内部パラメータを推定する、推定ステップと、
を含むことを特徴とする。

本発明の推定方法では、前記推定ステップにおいて、前記軌跡画像においてそれぞれ放物線状に現れる前記各軌跡の頂点を求め、前記各軌跡の頂点に基づいて、前記内部パラメータにおける、前記魚眼カメラの魚眼画像の歪みの中心を表す画像中心の座標を、推定すると、好適である。

本発明の推定方法では、前記推定ステップにおいて、前記軌跡画像中の前記各観測点にそれぞれ対応する前記特徴点の方位角及び仰角を、前記軌跡画像に基づいて求め、求めた前記各特徴点の前記方位角及び仰角を用いて、前記内部パラメータにおける、前記魚眼カメラの魚眼画像の歪みを表す歪みパラメータを、推定すると、好適である。

本発明の推定装置は、魚眼レンズを備えた魚眼カメラの内部パラメータの推定装置であって、
前記魚眼カメラが、複数の視覚的な特徴点を含む空間内で、該魚眼カメラの内部を通るとともに該魚眼カメラの光軸に対してほぼ垂直に延在する回転軸線の周りで回転される間に、該魚眼カメラによって複数の異なるタイミングのそれぞれで前記複数の特徴点を撮像することにより得られた、複数の魚眼画像を、取得する、取得部と、
前記各特徴点に対応する前記各魚眼画像中の観測点の集合から得られる軌跡を含む、軌跡画像を生成する、軌跡画像生成部と、
前記軌跡画像に基づいて、前記魚眼カメラの内部パラメータを推定する、推定処理部と、
を備えたことを特徴とする。

本発明の推定装置において、前記魚眼カメラを前記回転軸線の周りで回転させるように構成された、カメラ回転部をさらに備えていると、好適である。

本発明の推定プログラムは、魚眼レンズを備えた魚眼カメラの内部パラメータの推定プログラムであって、
前記魚眼カメラが、複数の視覚的な特徴点を含む空間内で、該魚眼カメラの内部を通るとともに該魚眼カメラの光軸に対してほぼ垂直に延在する回転軸線の周りで回転される間に、前記魚眼カメラによって複数の異なるタイミングのそれぞれで前記複数の特徴点を撮像することにより得られた、複数の魚眼画像を、取得する、取得ステップと、
前記各特徴点に対応する前記各魚眼画像中の観測点の集合から得られる軌跡を含む、軌跡画像を生成する、軌跡画像生成ステップと、
前記軌跡画像に基づいて、前記魚眼カメラの内部パラメータを推定する、推定ステップと、
をコンピュータに実行させることを特徴とする。

本発明によれば、簡単かつ高精度に、魚眼カメラの内部パラメータを推定できる、魚眼カメラの内部パラメータの、推定方法、推定装置、及び推定プログラムを提供することができる。

本発明の一実施形態に係る推定装置を示す概要図である。 魚眼カメラ座標系を説明するための図である。 魚眼画像座標系を説明するための図である。 本発明の一実施形態に係る推定方法を示すフローチャートである。 撮像ステップを説明するための図である。 軌跡画像の一例を示す図である。 推定ステップの一例を説明するためのフローチャートである。 シミュレーション結果を示す図である。

以下に図面を参照しつつ、本発明の実施形態について例示説明する。

（推定装置）
まず、図１を参照して、本発明の一実施形態に係る推定装置の一例を説明する。本実施形態の推定装置１は、魚眼レンズ１１を備えた魚眼カメラ１０の内部パラメータを推定するものである。図１の例において、推定装置１は、魚眼カメラ１０を支持するカメラ支持部９０と、コンピュータ７０と、を備えている。

本例において、魚眼カメラ１０は、ユーザによってシャッターが押された直後又は所定時間後に、撮影を開始し、撮影中は、複数の異なるタイミングで撮像して魚眼画像を生成し、ユーザによってシャッターが再び押された直後に、撮影を終了するように、設定されている。魚眼カメラ１０が撮影中に撮像するタイミングは、本例では、ある一定の撮像時間間隔毎（例えば0.033秒毎）に設定されている。
ただし、魚眼カメラ１０は、撮影を開始してから、所定の撮影時間（例えば10秒）の経過後に、自動的に撮影を終了するように、設定されてもよい。あるいは、魚眼カメラ１０の撮影の開始及び終了は、コンピュータ７０によって制御されてもよい。
また、魚眼カメラ１０は、ユーザによってシャッターが押される度に撮像するように構成されてもよい。

魚眼カメラ１０によって得られる時系列の魚眼画像は、後述するコンピュータ７０の記憶部７３に保存される。
より具体的には、本例において、魚眼カメラ１０は、撮影中又は撮影後に、無線通信又は有線通信によって、撮像により得られた時系列の魚眼画像を、後述するコンピュータ７０の通信部７１へ送信する。一方、コンピュータ７０側では、通信部７１が魚眼画像を受信すると、その魚眼画像が記憶部７３に格納される。
ただし、魚眼カメラ１０が、撮像により得られた時系列の魚眼画像を、魚眼カメラ１０に接続された外部メモリ（ＳＤカードメモリやＵＳＢメモリ等）に保存し、その後、ユーザが該外部メモリをコンピュータ７０に接続して、該外部メモリに保存された時系列の魚眼画像を、後述するコンピュータ７０の記憶部７３内に取り込むようにしてもよい。

本例において、カメラ支持部９０は、所定高さにある水平面内で、魚眼カメラ１０を回転させられるように、構成されている。より具体的に、カメラ支持部９０は、魚眼カメラ１０を上に載せるための回転台２０と、回転台２０の下端に固定され、鉛直方向に延びる回転軸３０と、回転軸３０の下部に出力軸が取り付けられ、回転軸３０を回転させるように構成された、モータ４０と、モータ４０に内蔵され、モータ４０の出力軸の回転角度位置（ひいては回転軸３０の回転角度位置）を出力する、エンコーダ５０と、モータ４０の下部に設けられた三脚６０と、を有している。
モータ４０の出力軸、回転軸３０及び回転台２０は、連動する。すなわち、モータ４０の出力軸の回転によって、回転軸３０及び回転台２０が回転され、これにより、回転台２０の上に設置された魚眼カメラ１０が回転される。

本例において、回転台２０の上面は水平方向と平行にされており、回転台２０の上面に載せられた魚眼カメラ１０の光軸Ｏも、水平方向と平行に向けられている。すなわち、光軸Ｏは水平面内にある。また、回転軸３０の中心軸線である回転軸線ＲＡは、鉛直方向に延びている。
魚眼カメラ１０は、回転軸３０の回転軸線ＲＡが魚眼カメラ１０の内部を通るように、回転台２０の上に載せられている。ここで、図１に示すように、回転軸３０の回転軸線ＲＡが、魚眼カメラ１０の内部で、魚眼カメラ１０の光軸Ｏと交差（より具体的には、直交）することが好ましい。

本例において、モータ４０は、後述するコンピュータ７０の回転制御部７２１によって、駆動制御される。本例では、カメラ支持部９０と、回転制御部７２１とが、魚眼カメラ１０を回転軸線ＲＡの周りで回転させる、カメラ回転部を構成している。
ただし、モータ４０は、コンピュータ７０の回転制御部７２１と同様の機能を持つ回転制御部を、モータ４０の内部に有していてもよい。

図１の例において、コンピュータ７０は、通信部７１と、制御部７２と、記憶部７３と、入力部７４と、表示部７５と、を有する。
なお、図１では、コンピュータ７０を機能ブロックにより模式的に図示しているにすぎず、コンピュータ７０のハードウェア構成は、通信部７１、制御部７２、記憶部７３、入力部７４、表示部７５の各機能が実現される限り、任意のものを採用してよい。例えば、コンピュータ７０は、そのハードウェア構成を観た場合に、１つの装置から構成されてもよいし、複数の装置から構成されてもよい。

本例において、通信部７１は、魚眼カメラ１０、モータ４０、及びエンコーダ５０との間で通信を行う。
なお、本例では、通信部７１が、魚眼カメラ１０から複数の魚眼画像を受信（取得）する、取得部を構成する。ただし、コンピュータ７０が、上述のように外部メモリ（ＳＤカードやＵＳＢメモリ）を介して、魚眼カメラ１０からの複数の魚眼画像を取得する場合、外部メモリと接続されるインタフェース（図示せず）が、取得部を構成する。

制御部７２は、例えば１つ又は複数のＣＰＵから構成され、記憶部７３に記憶された各種のプログラムを実行することにより、通信部７１、記憶部７３、入力部７４、及び表示部７５を含む、コンピュータ７０の全体を制御しながら、様々な処理を実行する。
制御部７２は、回転制御部７２１と、軌跡画像生成部７２２と、推定処理部７２３と、を有する。

回転制御部７２１は、記憶部７３に記憶された回転制御プログラム７３１を実行しながら、通信部７１を介して得られるエンコーダ５０の出力に基づいて、モータ４０の回転（ひいては魚眼カメラ１０の回転）を制御する。本例において、回転制御部７２１は、魚眼カメラ１０が撮影を行う間、モータ４０を、一定の回転速度で所定の回転角度（例えば180°）分だけ回転させる。
軌跡画像生成部７２２は、記憶部７３に記憶された軌跡画像生成プログラム７３２を実行しながら、魚眼カメラ１０から得られる時系列の魚眼画像に基づいて、後述する軌跡画像を生成する。
推定処理部７２３は、記憶部７３に記憶された推定処理プログラム７３３を実行しながら、軌跡画像生成部７２２が生成する軌跡画像に基づいて、魚眼カメラ１０の内部パラメータを推定する。
制御部７２による処理の詳細については、後にさらに詳しく説明する。

記憶部７３は、例えば１つ又は複数のＲＯＭ及び／又はＲＡＭから構成され、制御部７２が実行するための各種のプログラム（回転制御プログラム７３１、軌跡画像生成プログラム７３２、推定処理プログラム７３３等）や、魚眼カメラ１０から得られる時系列の魚眼画像等、様々な情報を記憶する。
なお、軌跡画像生成プログラム７３２と推定処理プログラム７３３とが、本実施形態の推定プログラムを構成している。

入力部７４は、例えばキーボード、マウス、及び／又は押しボタン等から構成され、ユーザからの入力を受け付ける。
表示部７５は、例えば液晶パネル等から構成され、後述する軌跡画像等、様々な情報を表示する。
なお、入力部７４及び表示部７５は、タッチパネルを構成してもよい。

（内部パラメータ）
ここで、図２及び図３を参照して、本発明の一実施形態に係る推定装置、推定方法、及び推定プログラムにより求められる、魚眼カメラの内部パラメータの一例について説明する。
図２は、魚眼カメラ１０を含む実空間の座標系（以下、「魚眼カメラ座標系」という。）を示している。図２に示す魚眼カメラ座標系は、魚眼カメラ１０の光軸ＯをＹ軸とし、魚眼カメラ１０の内部を通りＹ軸に直交する平面をＸＺ平面とする、３次元座標系である。ＸＹ平面は、水平面に平行である。魚眼カメラ座標系において、実空間上の点Ｐの方位角、すなわち、点Ｐから魚眼カメラ１０への投射線の、ＸＹ平面上への投影線が、Ｙ軸に対してなす角度を、αとする。また、点Ｐの仰角、すなわち、点Ｐから魚眼カメラ１０への投射線の、ＸＹ平面上への投影線が、該投射線に対してなす角度を、βとする。
このとき、点Ｐから魚眼カメラ１０への投射線と魚眼カメラ１０の光軸Ｏとのなす角θは、
で表される。

図３は、魚眼カメラ１０の魚眼画像の座標系（以下、「魚眼画像座標系」という。）を示している。図３に示す魚眼画像座標系は、魚眼カメラ１０の光軸Ｏと魚眼画像との交点である画像中心Ｃを原点とする、ｘｙ直交座標系である。魚眼画像座標系のｘ軸は、魚眼カメラ座標系のＸＹ平面と平行であり、魚眼画像座標系のｙ軸は、魚眼カメラ座標系のＺ軸と平行である。

実空間上の点Ｐは、魚眼カメラ１０によって、魚眼画像座標系における点ｐ（ｒ，φ）に投影される。
魚眼画像座標系の原点（画像中心Ｃ）と投影点ｐとを結ぶ線の、ｙ軸とのなす角φは、
で表わされる。

本例では、魚眼画像座標系の原点（画像中心Ｃ）から投影点ｐまでの距離ｒを、
で表すものとする。
ｋ₁、ｋ₃、ｋ₅は、魚眼カメラ１０の魚眼画像の歪みを表す、歪みパラメータである。ｋ₁は、１次の歪みパラメータであり、いわゆる焦点距離に関係している。ｋ₃は、３次の歪みパラメータである。ｋ₅は、５次の歪みパラメータである。

そして、式（３）に基づき、魚眼カメラ１０の内部パラメータＩは、
で表わされる。
ｃ_u、ｃ_vは、図３に示すように、魚眼画像の左上端を原点とするｕｖ直交座標系における、画像中心Ｃの座標である。画像中心Ｃは、魚眼カメラ１０の魚眼画像の歪みの中心を表す。
このように、本例における内部パラメータＩは、歪みパラメータｋ₁、ｋ₃、ｋ₅と、画像中心Ｃの座標ｃ_u、ｃ_vと、を含むものである。

なお、投影点ｐのｕｖ直交座標系における位置ベクトルをｍ（ｕ，ｖ）とすると、
となる。
以上の式（１）〜（６）により、投影点ｐの位置ベクトルｍは、α、β、Ｉで表わすことができる。

なお、魚眼画像座標系の原点（画像中心Ｃ）から投影点ｐまでの距離ｒは、式（３）以外の式から求めるようにしてもよい。式（３）により表される距離ｒは、１次、３次、及び５次の項からなるが、例えば、これに加えて、２次、４次、７次、及び／又は９次の項を加えてもよい。あるいは、式（３）における３次及び／又は５次の項を省略してもよい。この場合、距離ｒの式に応じて、内部パラメータＩの式も変わることとなる。

（推定方法）
つぎに、図４〜図８を参照して、本発明の一実施形態に係る魚眼カメラの内部パラメータの推定方法を説明する。図４は、本発明の一実施形態に係る推定方法を示すフローチャートである。ここでは、図１の例の推定装置１を用いる場合について説明するが、後述するように、他の構成を持つ推定装置を用いて、本実施形態に係る推定方法を実施することもできる。また、以下に説明する例では、推定の対象とする内部パラメータＩとして、式（４）に表されるものを用いるが、式（４）以外の式により表される内部パラメータＩを、本実施形態に係る推定方法によって推定するようにしてもよい。また、以下に説明する例では、図２及び図３を参照して説明した魚眼カメラ座標系及び魚眼画像座標系を用いるものとする。

−撮像ステップ（Ｓ１）−
まず、予め、魚眼カメラ１０が、カメラ支持部９０の回転台２０の上に設置される。本例において、回転台２０の回転軸線ＲＡは、魚眼カメラ１０の内部を通るとともに鉛直方向に延びており、魚眼カメラ１０の内部で、魚眼カメラ１０の光軸Ｏと直交している。
そして、図５に示すように、魚眼カメラ１０が、複数（この例では、４つ）の視覚的な特徴点８１ａ〜８４ａを含む空間内で、回転軸線ＲＡの周りで回転される間に、魚眼カメラ１０によって、複数の異なるタイミングのそれぞれで（ひいては、魚眼カメラ１０が複数の異なる回転位置のそれぞれにあるときに）、複数の特徴点８１ａ〜８４ａを撮像することにより、該複数の異なるタイミングに対応した複数の魚眼画像を得る（撮像ステップＳ１）。図５（ａ）及び図５（ｂ）は、それぞれ、撮像ステップＳ１において魚眼カメラ１０が回転軸線ＲＡの周りで回転されながら撮影する様子を示す、上面図及び側面図である。

より具体的に、本例では、コンピュータ７０の回転制御部７２１が、回転制御プログラム７３１に従い、通信部７１を介して得られるエンコーダ５０の出力に基づいて、モータ４０の回転を制御することにより、魚眼カメラ１０を、回転軸線ＲＡの周りで、一定の回転速度で180°分だけ回転させる。このときの魚眼カメラ１０の回転速度や回転角度（180°）は、予めユーザ等により設定される。
本例において、魚眼カメラ１０は、水平面に平行な仮想回転平面ＲＰ内で回転される。
また、魚眼カメラ１０は、回転軸線ＲＡの周りで回転される間、ある一定の撮像時間間隔毎（例えば0.033秒毎）に、その周囲にある特徴点８１ａ〜８４ａを撮像する。
魚眼カメラ１０は、撮影中又は撮影後に、無線通信又は有線通信によって、撮像により得られた時系列の魚眼画像を、コンピュータ７０の通信部７１へ送信する。一方、コンピュータ７０側では、通信部７１が魚眼画像を受信すると、制御部７２が、魚眼画像を記憶部７３に格納する。

視覚的な特徴点８１ａ〜８４ａとは、それが魚眼カメラ１０によって撮像されて魚眼画像に写ったときに、視覚的に識別し得るような点であり、例えば、魚眼カメラ１０の周りにある物体８１〜８４の角である。通常の環境内には、このような視覚的な特徴点が多数含まれているので、撮像ステップＳ１を実施するにあたっては、特別な被写体を用意する必要は無い。なお、図５に示す例では、４つの特徴点８１ａ〜８４ａを示しているが、特徴点は２個以上あればよい。また、特徴点の数が多いほど、より高精度に魚眼カメラの内部パラメータを推定できるので、好ましい。

なお、本例に限られず、回転軸線ＲＡは、魚眼カメラ１０の内部を通るとともに魚眼カメラ１０の光軸Ｏに対してほぼ垂直に延在していればよく、例えば鉛直方向と交わる方向に延びていてもよい。回転軸線ＲＡは、魚眼レンズ１１の内部（特に好ましくは、魚眼レンズ１１の中心（主点））を通るとともに魚眼カメラ１０の光軸Ｏに対してほぼ垂直に延在していることが、好ましい。なお、回転軸線ＲＡは、特徴点までの距離に比べて、魚眼レンズ１１の主点までの距離が、小さいほど良い。
また、魚眼カメラ１０による撮影と並行して行われる魚眼カメラ１０の回転の回転角度は、本例のように180°である必要はなく、180°より小さくてもよいし、180°より大きくてもよい。また、該回転角度は、180°以上であることが好ましく、360°であることがより好ましい。

−軌跡画像生成ステップ（Ｓ２）−
撮像ステップＳ１の後、図６に示すように、各特徴点８１ａ〜８４ａに対応する各魚眼画像中の観測点（像点）の集合から得られる軌跡Ｔ１〜Ｔ４を含む、軌跡画像を生成する（軌跡画像生成ステップＳ２）。図６は、軌跡画像の一例を模式的に表す図である。
より具体的には、まず、コンピュータ７０の軌跡画像生成部７２２が、軌跡画像生成プログラム７３２に従って、記憶部７３から時系列の魚眼画像を読み出す（取得する）。そして、軌跡画像生成部７２２は、軌跡画像生成プログラム７３２に従って、時系列の魚眼画像に基づいて、軌跡画像を生成する。ここで、例えば、軌跡画像生成部７２２は、各魚眼画像どうしを重ね合わせることにより、軌跡画像を生成する。あるいは、軌跡画像生成部７２２は、各魚眼画像から、それぞれの魚眼画像内に現れる特徴点８１ａ〜８４ａの観測点を抽出し、抽出した各観測点を１つの魚眼画像上に重ね合わせることにより、軌跡画像を生成する。その後、軌跡画像生成部７２２は、生成した軌跡画像を、記憶部７３に記憶させる。

なお、「観測点の集合から得られる軌跡」とは、観測点の集合そのものからなる軌跡と、観測点の集合に対して画像処理を施すことにより得られる、例えば、各観測点どうしを線分で結び、それらを連ねた線、あるいは、全観測点を滑らかに連ねた線（近似線）で示される軌跡との、いずれでもよい。本例では、観測点の集合そのものを軌跡としている。但し、図６では、説明の便宜のため、軌跡Ｔ１〜Ｔ４上の各観測点を滑らかに連ねた線を、破線で示している。

軌跡画像に現れる軌跡Ｔ１〜Ｔ４は、それぞれ別々の特徴点８１ａ〜８４ａに対応する。図６では、説明の便宜のため、各軌跡Ｔ１〜Ｔ４上の各観測点のｕｖ座標系での位置ベクトルを、ｍ_f1、ｍ_f2、ｍ_f3・・・ｍ_fN-2、ｍ_fN-1、ｍ_fNで表している。Ｎは、各軌跡Ｔ１〜Ｔ４上の観測点の数であり、撮像ステップＳ１で魚眼カメラ１０が撮像した回数（ひいては、魚眼画像の数）に等しい。

−推定ステップ（Ｓ３）−
軌跡画像生成ステップＳ２の後、軌跡画像生成ステップＳ２で得られた軌跡画像に基づいて、魚眼カメラ１０の内部パラメータを推定する（推定ステップＳ３）。
より具体的に、本例では、コンピュータ７０の推定処理部７２３が、推定処理プログラム７３３に従い、軌跡画像に基づいて、式（４）で表される内部パラメータＩを推定する。
図７は、推定ステップＳ３において、推定処理部７２３が行う処理を説明するためのフローチャートである。
図７に示すように、推定処理部７２３は、内部パラメータＩを推定するにあたって、まず、内部パラメータＩのうちの画像中心Ｃの座標ｃ_u、ｃ_vを、軌跡画像を用いて推定する（ステップＳ３１）。
つぎに、実際の軌跡画像中の各観測点に対応する実空間上の点の方位角α、仰角βの初期値を、軌跡画像を用いて、大まかに推定する（ステップＳ３２）。その後は、内部パラメータＩのうちの歪みパラメータｋ₁、ｋ₃、ｋ₅を、軌跡画像から求まる方位角α、仰角βを用いて、推定する（ステップＳ３３〜Ｓ５５）。ここで、推定処理部７２３は、いったん推定された歪みパラメータｋ₁、ｋ₃、ｋ₅を用いたシミュレーションによって、実際の軌跡画像中の観測点に対応する実空間上の点を魚眼画像へ投影した場合に得られる点（以下、「再投影点」という。）が、実際の観測点になるべく近づくように、修正パウエル法等の最適化手法により、歪みパラメータｋ₁、ｋ₃、ｋ₅の推定を行う。

以下、図７の各ステップＳ３１〜Ｓ５５での処理について、順番に、詳しく説明する。
まず、推定処理部７２３は、画像中心Ｃの座標ｃ_u、ｃ_vを推定する（ステップＳ３１）。軌跡画像中の軌跡Ｔ１〜Ｔ４は、画像中心Ｃを通るとともにｖ軸に平行な線（ｕ＝ｃ_uの線）に対して略線対称な放物線状に現れるという特徴がある。また、軌跡画像中の軌跡Ｔ１〜Ｔ４は、画像中心Ｃに近い軌跡ほど、軌跡の形状が緩やかである（傾きが小さい）という特徴がある。画像中心Ｃの座標ｃ_u、ｃ_vを推定するにあたっては、このような軌跡画像の特徴を利用する。
より具体的に、推定処理部７２３は、軌跡画像（図６）においてそれぞれ放物線状に現れる各軌跡Ｔ１〜Ｔ４の頂点を求め、各軌跡Ｔ１〜Ｔ４の頂点に基づいて、画像中心Ｃの座標ｃ_u、ｃ_vを推定する。例えば、軌跡画像における各軌跡Ｔ１〜Ｔ４を、図６において破線で示されたような放物線（二次曲線）で、近似する。つぎに、放物線状に現れる各軌跡Ｔ１〜Ｔ４の頂点のｕ座標を計算し、その平均を、画像中心Ｃのｕ座標ｃ_uとする。また、各軌跡Ｔ１〜Ｔ４のうち、近似された放物線の形状が最も緩やかである（放物線の傾きが最も小さい）ものを選択し、選択した軌跡（図６の例では、軌跡Ｔ１）の頂点のｖ座標を、画像中心Ｃのｖ座標ｃ_vとする。

つぎに、推定処理部７２３は、軌跡画像中の各観測点に対応する実空間上の点の方位角α、仰角βの初期値を、大まかに推定する（ステップＳ３２）。
まず、軌跡画像中の軌跡Ｔ１〜Ｔ４毎に、軌跡画像における画像中心Ｃを通るとともにｖ軸に平行な線（ｕ＝ｃ_uの線）上又はそれの最も近くにある観測点を特定し、特定した各観測点にそれぞれ対応する仰角βを計算する。ここで、ｕ＝ｃ_uの線上又はそれの最も近くにある各観測点が等距離射影されたものであるとの仮定のもと、画像中心Ｃからの距離に比例して仰角βの絶対値が大きくなり、軌跡画像の上端から下端までの間で仰角βが180°分だけ変化するものとして、仰角βを計算する。また、ｕ＝ｃ_uの線上又はそれの最も近くにある各観測点に対応する方位角αを、α＝０とする。このようにして、各軌跡Ｔ１〜Ｔ４におけるｕ＝ｃ_uの線上又はそれの最も近くにある観測点について、それぞれ方位角α及び仰角βの初期値を、大まかに推定する。
つぎに、ｕ＝ｃ_uの線上又はそれの最も近くにある観測点以外の各観測点に対応する方位角αを計算する。ここで、軌跡Ｔ１〜Ｔ４毎に観たときに、１つの軌跡の左端から右端に向かう間、方位角αの値は、軌跡の左端から軌跡に沿って測った距離に比例して変化し、軌跡の右端に至るまでには、撮像ステップＳ１で魚眼カメラ１０が回転された回転角度である180°分だけ変化するものと、仮定する。また、１つの軌跡上の各観測点にそれぞれ対応する仰角βは、同じである。よって、ｕ＝ｃ_uの線上又はそれの最も近くにある観測点以外の各観測点にそれぞれ対応する仰角βは、ｕ＝ｃ_uの線上又はそれの最も近くにある該観測点について求めた仰角βと、同じである。このようにして、ｕ＝ｃ_uの線上又はそれの最も近くにある観測点以外の各観測点についても、方位角α及び仰角βの初期値を大まかに推定する。

その後、推定処理部７２３は、上述のように求めた、軌跡画像中の各観測点に対応する方位角α及び仰角βの初期値を用いて、以下の式（７）〜（９）で求められる評価関数Ｄを計算する。
ここで、ｍ_fiは、観測点のｕｖ座標系での位置ベクトルであり、ｍ_riは、再投影点のｕｖ座標系での位置ベクトルである。再投影点ｍ_riは、式（１）〜（６）により、α、β、Ｉの関数で表すことができる。ｉ＝１は、各軌跡の左端に対応する。Ｎは、１つの軌跡上の観測点の数である。Ｍは、魚眼画像内の軌跡の数（ひいては特徴点の数）であり、図６の例では、Ｍ＝４である。
評価関数Ｄは、観測点からなる軌跡と、再投影点からなる軌跡との一致度合いを表すものであり、Ｄ_rjとＤ_sjとの和である。Ｄ_rjは、再投影誤差、すなわち、１つの軌跡における、観測点と、その観測点に対応する再投影点との距離の総和を表す。Ｄ_sjは、観測点からなる軌跡と、再投影点からなる軌跡との、形状の差を表す値である。
ステップＳ３２において、推定処理部７２３は、ｋ₃＝０、ｋ₅＝０とおくとともに、ｋ₁の値を変えながら評価関数Ｄを計算していき、評価関数Ｄの最小値Ｄ_minを求める。

つぎに、推定処理部７２３は、歪みパラメータｋ₁、ｋ₃、ｋ₅のうち、ｋ₁を推定する（ステップＳ３３）。より具体的には、直前に求めたＤ_minのときのｋ₁を、その推定値とする。

つぎに、推定処理部７２３は、軌跡画像中の各観測点に対応する方位角α及び仰角βを、本例ではステップＳ３２とは異なる方法で、再び計算する（ステップＳ３４）。
まず、軌跡画像中の軌跡Ｔ１〜Ｔ４毎に、ｕ＝ｃ_uの線上又はそれの最も近くにある観測点を特定し、特定した各観測点にそれぞれ対応する仰角βを計算する。ここで、特定した各観測点の方位角αを、α＝０とおくと、式（１）から、θ＝βが得られる。このθ＝βを式（３）に代入する。また、式（３）において、ｋ₁はステップＳ３３での推定値を用い、ｋ₃＝０、ｋ₅＝０とおくとともに、ｒは軌跡画像から求めたものを用いる。すると、βの方程式が得られ、βが求まる。このようにして、各軌跡に対応するβが得られる。
つぎに、ｕ＝ｃ_uの線上又はそれの最も近くにある観測点以外の各観測点に対応する方位角αを、上記のとおり既知のｋ₁、ｋ₃、ｋ₅、βと軌跡画像から求めた該各観測点のｒとを用いて、式（１）、（２）によりθを消してαの方程式を得ることにより、求める。
このようにして、軌跡画像中の各観測点に対応する方位角α及び仰角βが求まる。
その後、推定処理部７２３は、上述のように再計算した方位角α及び仰角βを用いて、上記のＤ_minを再び計算する。

つぎに、推定処理部７２３は、今回算出したＤ_minが、前回算出したＤ_min以下であれば（ステップＳ３５、ＮＯ）、ステップＳ３３に戻って再びｋ₁の推定を行う。この場合、ステップＳ３３では、直前に求めたＤ_min（すなわちステップＳ３４で求めたＤ_min）のときのｋ₁を、その推定値とする。
以下、ステップＳ３２〜Ｓ３５の処理を、「ループ処理１」という。
一方、今回算出したＤ_minが前回算出したＤ_minよりも大きければ（ステップＳ３５、ＹＥＳ）、ループ処理１で求めたＤ_minのうち最小のＤ_min（即ち、前回算出したＤ_min）と、該最小のＤ_minのときの歪みパラメータｋ₁とを、いったん記憶部７３に記憶させる（ステップＳ３６）。

つぎに、推定処理部７２３は、歪みパラメータｋ₁、ｋ₃、ｋ₅のうち、ｋ₁、ｋ₃を推定する（ステップＳ３７）。ここで、初めてステップＳ３７を実行する場合、推定処理部７２３は、ｋ₅＝０とおくとともに、ｋ₁、ｋ₃の値を変えながら評価関数Ｄを計算していき、評価関数Ｄの最小値Ｄ_minを求める。ここで、方位角α及び仰角βは、直前にステップＳ３４で求めた値を用いる。そして、得られたＤ_minのときのｋ₁、ｋ₃の値を、それらの推定値とする。
なお、２回目以降にステップＳ３７を実行する場合は、直前に求めたＤ_minのときのｋ₁、ｋ₃を、それらの推定値とする。

つぎに、推定処理部７２３は、軌跡画像中の各観測点に対応する方位角α及び仰角βを、ステップＳ３４と同様の方法で、再計算する（ステップＳ３８）。ただし、ステップＳ３４とは異なり、ｋ₅＝０とおく。
その後、推定処理部７２３は、上記再計算した方位角α及び仰角βを用いて、上記Ｄ_minを再び計算する。

つぎに、推定処理部７２３は、ステップＳ３７、Ｓ３８の処理を１０回繰り返し、その後（ステップＳ３９、ＹＥＳ）、その１０回繰り返した間に、少なくとも１回、Ｄ_minが小さくならなかった場合（ステップＳ４０、ＹＥＳ）、後述のステップＳ４４に進む。
一方、上記１０回繰り返した間に、少なくとも１回、Ｄ_minが小さくなった場合（ステップＳ４０、ＹＥＳ）、さらにステップＳ３７、Ｓ３８と同様の処理を行う（ステップＳ４１、Ｓ４２）。ここで、ステップＳ４１を実行する場合は、上記１０回繰り返した間における最小のＤ_minのときのｋ₁、ｋ₃を、それらの推定値とする。ステップＳ４２の後、今回算出したＤ_minが、前回算出したＤ_min以下であれば（ステップＳ４３、ＮＯ）、再びステップＳ４１、Ｓ４２を行う。一方、今回算出したＤ_minが、前回算出したＤ_minよりも大きければ（ステップＳ４３、ＹＥＳ）、後述のステップＳ４４に進む。
以下、ステップＳ３７〜Ｓ４３の処理を、「ループ処理２」という。
ステップＳ４４では、ループ処理２で求めたＤ_minのうち最小のＤ_minと、該最小のＤ_minのときの歪みパラメータｋ₁、ｋ₃とを、いったん記憶部７３に記憶させる。

その後、ステップＳ４５において、直前の３回にわたって記憶された最小のＤ_min（すなわち直前のステップＳ４４、直前のステップＳ５３、及び１回前のステップＳ４４のそれぞれで記憶された最小のＤ_min）どうしを比較する。ステップＳ５３の処理については、後述する。そして、最小のＤ_minが、直前の２回（すなわち直前のステップＳ４４及び直前のステップＳ５３）で連続して大きくならなかった場合（ステップＳ４５、ＮＯ）、ステップＳ４６に進み、ループ処理２と同様の処理（以下、「ループ処理３」という。）を、歪みパラメータｋ₁、ｋ₃、ｋ₅の全てについて行う（ステップＳ４６〜Ｓ５２）。その後、ステップＳ５３において、ループ処理３で求めたＤ_minのうち最小のＤ_minと、該最小のＤ_minのときの歪みパラメータｋ₁、ｋ₃、ｋ₅とを、いったん記憶部７３に記憶させる。
一方、ステップＳ４５において、直前の２回で連続して大きくなった場合は（ステップＳ４５、ＹＥＳ）、後述のステップＳ５５に進む。

ステップＳ５３の後のステップＳ５４において、最小のＤ_minが、直前の２回（すなわち直前のステップＳ５３及び直前のステップＳ４４）で連続して大きくならなかった場合（ステップＳ５４、ＮＯ）、再びステップＳ３７に戻って、ループ処理２から処理を再び行う。一方、ステップＳ５４において、直前の２回で連続して大きくなった場合は（ステップＳ５４、ＹＥＳ）、ステップＳ５５に進む。
ステップＳ５５において、推定処理部７２３は、記憶部７３を参照して、それまでに求めたＤ_minのうち最小のＤ_minと、該最小のＤ_minのときの歪みパラメータｋ₁、ｋ₃、ｋ₅とを、最終的な推定結果として出力し、処理を終了する。

（効果）
図８は、上述した例に係る撮像ステップＳ１、軌跡画像生成ステップＳ２を経て生成された、観測点からなる軌跡Ｔ_fを含む軌跡画像に、上記の例における推定ステップＳ３におけるシミュレーションで得られた、再投影点からなる軌跡Ｔ_rを、重ね合わせて、軌跡画像の右上部分を拡大した図である。図８（ａ）は、推定ステップＳ３の実行途中に得られた歪みパラメータの推定値を含む内部パラメータＩを用いて、再投影点の軌跡Ｔ_rを得た場合の図であり、再投影点の軌跡Ｔｒが観測点の軌跡Ｔ_fから少しずれている。図８（ｂ）は、推定ステップＳ３の最終結果として得られた歪みパラメータの推定値を含む内部パラメータＩを用いて、再投影点の軌跡Ｔｒを得た場合の図であり、再投影点の軌跡Ｔ_rが観測点の軌跡Ｔ_fにほぼ重なっている。図８の結果から、上述した例に係る撮像ステップＳ１、軌跡画像生成ステップＳ２を用いて生成された軌跡画像を用いることにより、高い精度で内部パラメータＩを推定できることが判る。

本実施形態によれば、特別な被写体を用意することを要さず、簡単かつ高精度に、魚眼カメラの内部パラメータを推定することができる。

なお、本施形態の推定装置や本実施形態の推定方法（ひいては本実施形態の推定プログラム）は、それぞれ、様々な変更例が可能である。
例えば、撮像ステップＳ１において、ユーザが、魚眼カメラ１０の回転や、魚眼カメラ１０への撮像の指示（シャッターボタンを押す等）を、行ってもよい。
また、撮像ステップＳ１において、撮像時間間隔は、一定である必要はなく、まだらな間隔でもよい。

推定ステップＳ３では、上述した方法以外の方法によっても、軌跡画像に基づいて、内部パラメータＩを推定してもよい。
より具体的には、ステップＳ３１において、画像中心Ｃのｖ座標ｃ_vを求めるにあたっては、軌跡画像内の各軌跡のうち、近似された放物線の形状が最もフラット（すなわち最も直線に近い）ものを特定し、特定した軌跡上の各観測点のｖ座標の平均値を、画像中心Ｃのｖ座標ｃ_vとしてもよい。あるいは、軌跡画像内の各軌跡から近似された各放物線の頂点のｖ座標及び２次の項の係数を求め、２次の項の係数が０になるｖ座標を最小２乗法で推定して、その結果を、画像中心Ｃのｖ座標ｃ_vとしてもよい。

また、上記の例では、歪みパラメータの推定において、軌跡画像における観測点と再投影点とを１対１の対応関係とし、方位角α及び仰角βの算出や、評価関数Ｄの算出を行うものであるが、これに限られない。例えば、評価関数Ｄを、軌跡画像における観測点と、再投影点が描く軌跡との、距離の式で表し、方位角αを算出せずに、仰角βのみを算出するようにしてもよい。
また、歪みパラメータの推定において、軌跡画像における観測点を用いる代わりに、観測点が描く軌跡あるいは該軌跡から近似された放物線上の、任意の点を用いてもよい。

また、撮像ステップＳ１での各撮像タイミングにおける魚眼カメラ１０の回転位置から、方位角αを求めるようにしてもよい。その場合、魚眼カメラ１０の回転位置は、エンコーダ５０の出力に基づいて求められる。あるいは、モータ４０をステッピングモータにより構成することで、各撮像タイミングにおける魚眼カメラ１０の回転位置を予め設定したとおりに制御できるようにしてもよい。その場合、エンコーダ５０を省略してもよい。

本発明の推定方法、推定装置、及び推定プログラムは、魚眼カメラの内部パラメータを推定する際に利用できる。

１ 推定装置
１０ 魚眼カメラ
１１ 魚眼レンズ
２０ 回転台
３０ 回転軸
４０ モータ
５０ エンコーダ
６０ 三脚
７０ コンピュータ
７１ 通信部（取得部）
７２ 制御部
７３ 記憶部
７４ 入力部
７５ 表示部
８１〜８４ 物体
８１ａ〜８４ａ 特徴点
９０ カメラ支持部（カメラ回転部）
７２１ 回転制御部（カメラ回転部）
７２２ 軌跡画像生成部
７２３ 推定処理部
７３１ 回転制御プログラム
７３２ 軌跡画像生成プログラム（推定プログラム）
７３３ 推定処理プログラム（推定プログラム）
Ｃ 画像中心
Ｏ 光軸
ＲＡ 回転軸線
ＲＰ 仮想回転平面
Ｔ１〜Ｔ４ 軌跡

Claims (6)

1. 魚眼レンズを備えた魚眼カメラの内部パラメータの推定方法であって、
   前記魚眼カメラが、複数の視覚的な特徴点を含む空間内で、該魚眼カメラの内部を通るとともに該魚眼カメラの光軸に対してほぼ垂直に延在する回転軸線の周りで回転される間に、前記魚眼カメラによって、複数の異なるタイミングのそれぞれで前記複数の特徴点を撮像することにより、複数の魚眼画像を得る、撮像ステップと、
   前記各特徴点に対応する前記各魚眼画像中の観測点の集合から得られる軌跡を含む、軌跡画像を生成する、軌跡画像生成ステップと、
   前記軌跡画像に基づいて、前記魚眼カメラの内部パラメータを推定する、推定ステップと、
   を含むことを特徴とする、推定方法。
2. 前記推定ステップにおいて、前記軌跡画像においてそれぞれ放物線状に現れる前記各軌跡の頂点を求め、前記各軌跡の頂点に基づいて、前記内部パラメータにおける、前記魚眼カメラの魚眼画像の歪みの中心を表す画像中心の座標を、推定する、請求項１に記載の推定方法。
3. 前記推定ステップにおいて、前記軌跡画像中の前記各観測点にそれぞれ対応する前記特徴点の方位角及び仰角を、前記軌跡画像に基づいて求め、求めた前記各特徴点の前記方位角及び仰角を用いて、前記内部パラメータにおける、前記魚眼カメラの魚眼画像の歪みを表す歪みパラメータを、推定する、請求項１又は２に記載の推定方法。
4. 魚眼レンズを備えた魚眼カメラの内部パラメータの推定装置であって、
   前記魚眼カメラが、複数の視覚的な特徴点を含む空間内で、該魚眼カメラの内部を通るとともに該魚眼カメラの光軸に対してほぼ垂直に延在する回転軸線の周りで回転される間に、該魚眼カメラによって複数の異なるタイミングのそれぞれで前記複数の特徴点を撮像することにより得られた、複数の魚眼画像を、取得する、取得部と、
   前記各特徴点に対応する前記各魚眼画像中の観測点の集合から得られる軌跡を含む、軌跡画像を生成する、軌跡画像生成部と、
   前記軌跡画像に基づいて、前記魚眼カメラの内部パラメータを推定する、推定処理部と、
   を備えたことを特徴とする、推定装置。
5. 前記魚眼カメラを前記回転軸線の周りで回転させるように構成された、カメラ回転部をさらに備えた、請求項４に記載の推定装置。
6. 魚眼レンズを備えた魚眼カメラの内部パラメータの推定プログラムであって、
   前記魚眼カメラが、複数の視覚的な特徴点を含む空間内で、該魚眼カメラの内部を通るとともに該魚眼カメラの光軸に対してほぼ垂直に延在する回転軸線の周りで回転される間に、前記魚眼カメラによって複数の異なるタイミングのそれぞれで前記複数の特徴点を撮像することにより得られた、複数の魚眼画像を、取得する、取得ステップと、
   前記各特徴点に対応する前記各魚眼画像中の観測点の集合から得られる軌跡を含む、軌跡画像を生成する、軌跡画像生成ステップと、
   前記軌跡画像に基づいて、前記魚眼カメラの内部パラメータを推定する、推定ステップと、
   をコンピュータに実行させることを特徴とする、推定プログラム。