WO2023084601A1

WO2023084601A1 - 算出装置、算出方法、カメラパラメータ推定装置、及び非一時的なコンピュータ可読媒体

Info

Publication number: WO2023084601A1
Application number: PCT/JP2021/041187
Authority: WO
Inventors: 学中野
Original assignee: 日本電気株式会社
Priority date: 2021-11-09
Filing date: 2021-11-09
Publication date: 2023-05-19

Abstract

算出装置において算出部は、ラグランジュ乗数を用いずにＥ行列（Essential matrix）の局所解を表す、「所定の連立方程式」を解くことによって、Ｅ行列を算出する。この「所定の連立方程式」は、Ｅ行列の特徴を制約条件としたラグランジュ関数ＬであってＮ組の対応点ペアについてのエピポーラ誤差に対するラグランジュ関数Ｌの１次最適性条件を表している。Ｎは、例えば、６以上の自然数である。各「対応点ペア」は、第１画像及び第２画像にそれぞれ含まれ且つ互いに対応する２つの対応点を含む。

Description

算出装置、算出方法、カメラパラメータ推定装置、及び非一時的なコンピュータ可読媒体

　本開示は、算出装置、算出方法、及び非一時的なコンピュータ可読媒体に関する。

　複数台のカメラを用いて同一の物体（被写体）を撮影した複数の画像を含む画像列から、カメラパラメータとその物体の３次元情報とを復元する問題は、Structure-from-Motion（SfM）や多視点幾何学問題と呼ばれている。カメラパラメータは、「内部パラメータ」及び「外部パラメータ」の２種類を含む。内部パラメータは、焦点距離、レンズ歪、光学中心座標などのレンズ固有のパラメータであり、外部パラメータは、カメラ間の３次元的な回転行列及び並進ベクトルである。内部パラメータは、レンズが固定ならば事前に計測が可能であり、内部パラメータが既知のカメラは、校正済みカメラと呼ばれる。
　校正済みカメラによる２視点幾何において「カメラ間の相対的な外部パラメータ」は、３×３の基本行列（Ｅ行列：Essential matrix）として表現される。ここで「相対的」としているのは、絶対的な世界座標系の定義が不明な場合、１つ目のカメラの座標系を基準として２つ目のカメラの回転行列及び並進ベクトルを表現するためである。この相対的な並進ベクトルと回転行列とを乗算したものがＥ行列である。このＥ行列の特徴として、固有値の１つがゼロであり、２つの特異値が等しい、ことが知られている。

　Ｅ行列の算出方法として、D. Nisterの「５点法」が広く知られている。D. Nisterの「５点法」は、２枚の画像における２つの対応点を含む「対応点ペア」を、少なくとも５組用いる。各対応点ペアの２つの対応点は、同じ３次元座標がそれぞれのカメラで観測された２つの画像座標である。

　５点法では、最大で１０個の実数解が算出されるため、一意に決めることが困難なことがある。解を一意に定めるためには、６組以上の対応点ペアを用いて、最小自乗解を算出すればよい。その方法の一つとして、８組以上の対応点ペアを用いて一意の解を算出するR. Hartleyの８点法も広く用いられている。特許文献１には、８点法を用いて２視点幾何を多視点カメラに適用し、物体の３次元形状の復元を行う方法が開示されている。また、特許文献２には、５点法を用いて２視点幾何を多視点カメラに適用し、物体の３次元形状の復元を行う方法が開示されている。ここで、「２視点幾何」とは、SfMの最小条件である２つの校正済みカメラを用いることに対応する。

国際公開第２０１２／０１４４３０号国際公開第２０１９／０５８４８７号

　しかしながら、特許文献１及び特許文献２に開示の５点法及び８点法は、統計的に最適な一意のＥ行列を算出できない、可能性がある。

　本開示の目的は、統計的に最適な一意のＥ行列を算出できる、算出装置、算出方法、及び非一時的なコンピュータ可読媒体を提供することにある。

　第１の態様にかかる算出装置は、対応点ペアの２つの対応点間の幾何学的な制約であるエピポーラ制約を表すための基本行列（Ｅ行列）を算出する算出装置であって、
　複数の対応点ペアのそれぞれに含まれる各対応点の画像平面内座標を取得し、各対応点ペアが第１画像及び第２画像にそれぞれ含まれ且つ互いに対応する２つの対応点を含む、取得部と、
　前記Ｅ行列の特徴を制約条件としたラグランジュ関数Ｌであって前記複数の対応点ペアについてのエピポーラ誤差に対するラグランジュ関数Ｌの１次最適性条件を表し、且つ、ラグランジュ乗数を用いずに前記Ｅ行列の局所解を表す、所定の連立方程式を解くことによって、前記Ｅ行列を算出する算出部と、
　を具備する。

　第２の態様にかかる算出方法は、算出装置によって実行される、対応点ペアの２つの対応点間の幾何学的な制約であるエピポーラ制約を表すための基本行列（Ｅ行列）を算出する算出方法であって、
　複数の対応点ペアのそれぞれに含まれる各対応点の画像平面内座標を取得することを含み、
　各対応点ペアが第１画像及び第２画像にそれぞれ含まれ且つ互いに対応する２つの対応点を含み、
　前記算出方法は、前記Ｅ行列の特徴を制約条件としたラグランジュ関数Ｌであって前記複数の対応点ペアについてのエピポーラ誤差に対するラグランジュ関数Ｌの１次最適性条件を表し、且つ、ラグランジュ乗数を用いずに前記Ｅ行列の局所解を表す、所定の連立方程式を解くことによって、前記Ｅ行列を算出することを含む。

　第３の態様にかかるプログラムは、対応点ペアの２つの対応点間の幾何学的な制約であるエピポーラ制約を表すための基本行列（Ｅ行列）を算出する処理を、算出装置に実行させるプログラムであって、
　前記処理は、複数の対応点ペアのそれぞれに含まれる各対応点の画像平面内座標を取得することを含み、
　各対応点ペアが第１画像及び第２画像にそれぞれ含まれ且つ互いに対応する２つの対応点を含み、
　前記処理は、前記Ｅ行列の特徴を制約条件としたラグランジュ関数Ｌであって前記複数の対応点ペアについてのエピポーラ誤差に対するラグランジュ関数Ｌの１次最適性条件を表し、且つ、ラグランジュ乗数を用いずに前記Ｅ行列の局所解を表す、所定の連立方程式を解くことによって、前記Ｅ行列を算出することを含む。

　本開示により、統計的に最適な一意のＥ行列を算出できる、算出装置、算出方法、及び非一時的なコンピュータ可読媒体を提供することができる。

第１実施形態における算出装置の一例を示すブロック図である。２視点幾何の具体例の説明に供する図である。第２実施形態における算出装置の一例を示すブロック図である。第２実施形態における算出装置による処理動作の一例を示すフローチャートである。第３実施形態における算出装置の一例を示すブロック図である。第３実施形態における算出装置による処理動作の一例を示すフローチャートである。第４実施形態におけるカメラパラメータ推定装置の一例を示すブロック図である。算出装置のハードウェア構成例を示す図である。

　以下、図面を参照しつつ、実施形態について説明する。なお、実施形態において、同一又は同等の要素には、同一の符号を付し、重複する説明は省略される。

＜第１実施形態＞
　図１は、第１実施形態における算出装置の一例を示すブロック図である。図１に示す算出装置１０は、「対応点ペア」の２つの対応点間の幾何学的な制約である「エピポーラ制約」を表すための「基本行列（Ｅ行列：Essential matrix）」を算出する。

　図１において算出装置１０は、取得部１１と、算出部１２とを有している。

　取得部１１は、Ｎ組の「対応点ペア」のそれぞれに含まれる各対応点の画像平面内座標を取得する。画像平面内座標は、「観察画像座標」と呼んでもよい。Ｎは、例えば、６以上の自然数である。各「対応点ペア」は、例えば、「第１画像」及び「第２画像」にそれぞれ含まれ且つ互いに対応する２つの対応点を含む。

　ここで、「２視点幾何」の具体例について説明する。図２は、２視点幾何の具体例の説明に供する図である。図２には、カメラ１のカメラ座標系及びカメラ２のカメラ座標系が示されている。図２において画像ＰＣ１は、カメラ１によって撮影された画像であり、画像ＰＣ２は、カメラ２によって撮影された画像である。すなわち、例えば、画像ＰＣ１は、上記の「第１画像」に対応し、画像ＰＣ２は、上記の「第２画像」に対応する。

　そして、図２において、世界座標系における或る３次元点Ｘが、カメラ１によって画像平面内座標ｍとして観測されている。この画像平面内座標ｍは、例えば、画像座標ｍと呼んでもよいし、観察画像座標ｍと呼んでもよい。また、図２において、その３次元点Ｘが、カメラ２によって画像平面内座標ｍ’として観測されている。観察画像座標ｍ及び画像平面内座標ｍ’は、互いに対応する点であるので、それぞれ「対応点」である。そして、観察画像座標ｍ及び画像平面内座標ｍ’を纏めて「対応点ペア」と呼ぶことができる。

　本実施形態において世界座標系は明示的に与えられないため、３次元点Xの座標値は不明である。世界座標系は、任意に設定可能なため、カメラ１のカメラ座標系と一致しているものとする。すなわち、カメラ１の外部パラメータにおいて、位置（つまり、３次元並進ベクトル）は、[0,0,0]であり、回転行列は、３×３の単位行列である。そして、カメラ２の外部パラメータは、並進ベクトルt、回転行列Ｒで表される。なお、２つのカメラ１，２は校正済みであり、その内部パラメータは、既知である。

　取得部１１が対応点（対応点の画像平面内座標）を取得する方法として、例えばＳＩＦＴ（Scale Invariant Feature Transform）やＳＵＲＦ（Speeded Up Robust Features）といった広く知られている手法を利用してもよい。ＳＩＦＴやＳＵＲＦは誤った対応点を取得することがあるため、例えば広く知られている５点法とＲＡＮＳＡＣ（Random Sample Consensus）とを組み合わせて誤った対応点を除去する方法を用いてもよい。または、取得部１１は、ユーザの手動による、対応点の指定や誤対応点の除去を受け付けることによって、対応点を取得してもよい。

　算出部１２は、ラグランジュ乗数を用いずにＥ行列の局所解を表す、「所定の連立方程式」を解くことによって、Ｅ行列を算出する。この「所定の連立方程式」は、さらに、Ｅ行列の特徴を制約条件としたラグランジュ関数Ｌであって上記Ｎ組の対応点ペアについてのエピポーラ誤差に対するラグランジュ関数Ｌの１次最適性条件を表してもよい。「エピポーラ誤差」は、例えば、Ｅ行列と対応点ペアとで構成されるエピポーラ制約（エピポーラ方程式）の代数学的誤差である。

　以上のように第１実施形態によれば、算出装置１０において算出部１２は、ラグランジュ乗数を用いずにＥ行列の局所解を表す、「所定の連立方程式」を解くことによって、Ｅ行列を算出する。この「所定の連立方程式」は、Ｅ行列の特徴を制約条件としたラグランジュ関数Ｌであって上記Ｎ組の対応点ペアについてのエピポーラ誤差に対するラグランジュ関数Ｌの１次最適性条件を表してもよい。

　この算出装置１０の構成により、統計的に最適な一意のＥ行列を算出できる。すなわち、算出装置１０の構成により、Ｅ行列と対応点ペアとで構成されるエピポーラ制約（エピポーラ方程式）の代数学的誤差（つまり、「エピポーラ誤差」）の局所解を一意に算出できる。

　なお、上記の算出装置１０は、対応点ペアの２つの対応点間の幾何学的な制約であるエピポーラ制約を表すための基本行列（Ｅ行列）を算出する算出方法であって、複数の対応点ペアのそれぞれに含まれる各対応点の画像平面内座標を取得することを含み、各対応点ペアが第１画像及び第２画像にそれぞれ含まれ且つ互いに対応する２つの対応点を含み、その算出方法は、上記Ｅ行列の特徴を制約条件としたラグランジュ関数Ｌであって上記複数の対応点ペアについてのエピポーラ誤差に対するラグランジュ関数Ｌの１次最適性条件を表し、且つ、ラグランジュ乗数を用いずに上記Ｅ行列の局所解を表す、所定の連立方程式を解くことによって、上記Ｅ行列を算出することを含む、算出方法を実行している。

　＜変形例＞
　以上の説明では、「第１画像」と「第２画像」とが、ある時刻において同一物体を、異なる２つのカメラによって異なる視点から撮影した、２枚の静止画であるものとして説明を行ったが、本実施形態は、これに限定されない。例えば、「第１画像」と「第２画像」とは、時系列の連続した動画の2枚フレーム画像でもよい。このとき、「第１画像」を撮影するカメラと「第２画像」を撮影するカメラとは、同一であってもよいし、異なっていてもよい。要するに、「第１画像」と「第２画像」とは、内部パラメータが校正済みの１つ以上のカメラによって、同一の物体または同一のシーンが異なる視点から撮影された画像であればよい。

＜第２実施形態＞
　第２実施形態は、より具体的な実施形態に関する。

　＜算出装置の構成例＞
　図３は、第２実施形態における算出装置の一例を示すブロック図である。図３において算出装置２０は、取得部１１と、算出部２１とを有している。算出部２１は、初期値決定部２１Ａと、局所解算出部２１Ｂとを有している。

　初期値決定部２１Ａは、取得部１１からＮ組の「対応点ペア」を受け取る。そして、初期値決定部２１Ａは、受け取ったＮ組の「対応点ペア」を用いて、Ｅ行列の初期値を決定する。この初期値は、「対応点ペア」に対して５点法や８点法を適用することによって算出されてもよい。なお、初期値決定部２１Ａは、取得部１１で取得された「対応点ペア」を用いずに、単にランダムな値を用いてＥ行列の初期値を決定してもよい。また、初期値決定部２１Ａは、「第１画像」と「第２画像」とが連続する動画の２フレームでありそれより前のフレームのＥ行列が既知である場合、その既知のＥ行列を初期値としてもよい。

　局所解算出部２１Ｂは、上記のＥ行列の初期値を用いて、「所定の連立方程式」を満たすように、Ｅ行列の最適化を実行する。「所定の連立方程式を満たす最適化」としては、例えばガウスニュートン法といった様々な公知の方が用いられてもよい。ここで、「所定の連立方程式」は、上記Ｅ行列の特徴を制約条件としたラグランジュ関数Ｌであって上記複数の対応点ペアについてのエピポーラ誤差に対するラグランジュ関数Ｌの「１次最適性条件」を表し、且つ、ラグランジュ乗数を用いずに上記Ｅ行列の局所解を表す非線形連立方程式である。本実施形態において「１次最適性条件」は、例えば、行列［Ｍｅ］に対して左からＥ行列の転置行列を掛けることで得られる行列及び行列［Ｍｅ］に対して右からＥ行列の転置行列を掛けることで得られる行列が対称行列となる「第１条件」、及び、Ｅ行列の特徴を表す「第２条件」を含む。上記の行列［Ｍｅ］は、対応点ペアの第１画像平面内座標ｍ及び第２画像平面内座標ｍ’から算出される係数行列Ｍと上記Ｅ行列を変形したベクトルｅとの、積であるベクトルを変形した行列である。すなわち、本実施形態において「１次最適性条件」は、例えば、３×３の行列表現に変形されると対称行列となり、ラグランジュ乗数を含まずにＥ行列の特徴を表現する、ことを特徴とする。Ｅ行列の特徴は、固有値の１つがゼロであり、２つの特異値が等しい、ことである。

　このように、Ｅ行列の初期値を用いてＥ行列の特徴を表す非線形方程式を最適化することで、Ｅ行列の１つの局所解を算出できる。すなわち、統計的に最適な解であるＥ行列の計算が実現可能である。

　＜算出装置の動作例＞
　以上の構成を有する算出装置２０の処理動作の一例について説明する。図４は、第２実施形態における算出装置による処理動作の一例を示すフローチャートである。

　初期値決定部２１Ａは、Ｅ行列の初期値を決定する（ステップＳ１０１）。

　局所解算出部２１Ｂは、Ｅ行列の初期値を用いて、「所定の連立方程式」を満たすＥ行列を１つ算出する（ステップＳ１０２）。

　ここで、局所解算出部２１Ｂによる処理動作の具体例について説明する。以下では、図２を参照して図２に示した記号を用いて説明する。すなわち、世界座標系の原点は、カメラ１のカメラ座標系に一致させている。或る３次元点Ｘは、２つのカメラ画像においてそれぞれｍ、ｍ’として観測されている。２つのカメラ１，２は、校正済みであり、これらの内部パラメータは、既知である。このため、ｍ、ｍ’は、いわゆる正規化画像座標である。さらに、ｍ、ｍ’は、３次元の斉次座標表現として、つまり、２次元の画像座標に１を要素に加えた３次元ベクトルとして表される。なお、以下の説明において、上付きＴは、ベクトル又は行列の転置を表す。

　まず、解くべき問題の定義を行う。

　Ｅ行列は、カメラ１からカメラ２への相対的な並進ベクトルｔ＝［ｔｘ，ｔｙ，ｔｚ］^Ｔと、回転行列Ｒとを用いて、次の式（１）で表される。

　そして、Ｅ行列は、１つの固有値がゼロであり、２つの特異値が等しい、という特徴を持つ。この特徴は、次の式（２）で表されることが知られている。

　なお、tr()は、行列のトレースを表す。

　Ｅ行列及びi番目の対応点ペアは、次の式（３）で表されるエピポーラ方程式（エピポーラ制約）を満たすことが知られている。

　式（２）、（３）が示すように、Ｅ行列には、符号の不定性がある。すなわち、符号を反転した（－Ｅ）を式（２）、（３）に代入しても、成立する。符号の不定性を持つＥ行列は、２組のtと２組のＲの、計４組のｔ及びＲで表されることが知られている。

　６組以上の対応点ペアが与えられたとき、一般に対応点ペアはノイズを含むため、式（３）は厳密にはゼロにならない。それゆえ、最適なＥ行列は、次の式（４）の最適化問題の解となる。

　ここで、eは、Ｅ行列を９次元ベクトルへ変形した表現である。

は、自明解Ｅ＝０を避けるための制約条件である。また、Ｍは、対応点ペアｍ、ｍ’から算出される９×９の係数行列であり、次の式で計算される。

は、クロネッカ積である。式（４）の目的関数ｅ^ＴＭｅは、代数学的誤差（エピポーラ誤差）と呼ばれている。

　次に、上記の式（４）の局所解を計算する方法について説明する。

　式（４）に対するラグランジュ関数Ｌと、このラグランジュ関数ＬのＥ行列に対する導関数は、次の式（５）で表される。

　ここで、λは、スカラーのラグランジュ乗数であり、Sは、９つのラグランジュ乗数を成分とする３×３の行列であり、１／２は、便宜上の数字である。

　次に、ラグランジュ関数ＬのＥ行列に対する導関数に対して、左右からＥ^Ｔを掛けると、式（６）が得られる。

　ここで、［Ｍｅ］は、行列Ｍとベクトルeとの積である９次元ベクトルを３×３の行列に変形した表現である。

　式（６）より、Ｅ^Ｔ［Ｍｅ］及び［Ｍｅ］Ｅ^Ｔがいずれも対称行列となることが、それぞれの右辺から示されている。すなわち、Ａ＝Ｅ^Ｔ［Ｍｅ］及びＢ＝［Ｍｅ］Ｅ^Ｔが対称行列となり、且つ、式（４）の制約条件を満たす、Ｅ行列（つまり、「１次最適性条件」を満たすＥ行列）を求めればよい。これは、次の式（７）で表される、Ｅ行列に関する連立方程式として表される。すなわち、この連立方程式は、複数の対応点ペアについてのエピポーラ誤差に対するラグランジュ関数ＬのＥ行列に対する、導関数に基づいて導出されている。また、この連立方程式は、ラグランジュ乗数を用いずにＥ行列の局所解を表している。

　ここで、Ａ_ｉｊ、Ｂ_ｉｊは、それぞれ、行列Ａの（ｉ，ｊ）成分及び行列Ｂの（ｉ，ｊ）成分を表す。

　すなわち、ステップＳ１０２において、局所解算出部２１Ｂは、Ｅ行列の初期値を用いて、式（７）の「所定の連立方程式」を満たす解（Ｅ行列）を１つ算出してもよい。

　以上のように第２実施形態によれば、算出装置２０において算出部２１は、行列［Ｍｅ］に対して左からＥ行列の転置行列を掛けることで得られる行列と行列［Ｍｅ］に対して右からＥ行列の転置行列を掛けることで得られる行列とが対称行列となる「第１条件」、及び、Ｅ行列の特徴を表す「第２条件」を含む「１次最適性条件」に基づいて得られる、所定の連立方程式を解くことによって、Ｅ行列を算出する。行列［Ｍｅ］は、対応点ペアの第１画像平面内座標ｍ及び第２画像平面内座標ｍ’から算出される係数行列Ｍと前記Ｅ行列を変形したベクトルｅとの積であるベクトルを変形した行列である。

　この算出装置２０の構成により、統計的に最適なＥ行列を一意に算出可能である。その理由は、以下の通りである。例えば、式（７）は、ラグランジュ乗数を含まずにＥ行列のみで式（４）の１次最適性条件を表現した連立方程式である。そして、Ｅ行列の初期値を用いた最適化により算出した連立方程式の解は、式（４）の最適化問題の局所解となる。すなわち、６組以上の対応点ペアが与えられた場合に統計的に最適な解を一意に算出することが可能となる。

　＜比較例＞
　比較例として、５点法及び８点法に関して考察する。

　まず、５点法は、解が一意に定まらない。そのため、複数の対応点ペアから５組を選び、選んだ対応点ペアを用いてＥ行列を算出する、という操作を繰り返し行う必要がある。対応点は、一般に数百個から数千個は得られるため、すべての組み合わせを試すことは現実的ではない。さらに、このような方法を持ってしても、５点法はすべての対応点ペアを同時に用いて統計的に最適な解を一意に算出することはできない。ここで、「統計的に最適な解」とは、Ｅ行列と対応点ペアとで構成されるエピポーラ制約（エピポーラ方程式）の代数学的誤差（つまり、「エピポーラ誤差」）の局所解のことである。

　次に、５点法及び８点法では、対応点ペアの数が６組または７組の場合にＥ行列を算出できない。なぜならば、５点法は、５組のみを対象とした方法であり、８点法は、８組以上を対象とした方法のためである。

　次に、１つの平面で構成されたシーンに対しては、８点法によってはＥ行列を算出できない。なぜならば、８点法の算出に用いる９×９の係数行列がランク落ちするため、理論的に解が一意に定まらなくなるからである。そのため、８点法には、車載画像のように１つの平面（路面）が支配的なシーンに対しては、Ｅ行列の算出に失敗したり、Ｅ行列の推定精度が著しく低下するという課題がある。５点法は、１つの平面であってもＥ行列を算出可能であるが、上記の通り、解を一意に定めることができず、算出する解は統計的に最適ではない。

　次に、８点法で算出される一意の解は、統計的に最適ではない。８点法は、まず９×９の係数行列の最小固有値に対応する９次元の固有ベクトルを計算する。次に、固有ベクトルを３×３の行列に変形し、前述した２つの特徴（つまり、１つの固有値がゼロであり、２つの特異値が等しいこと）を満たすように、特異値分解を用いた補正を行う。この補正操作は、補正前後における３×３行列のフロベニウスノルムの差分が最小となることを目的としており、補正後のＥ行列は、統計的に最適ではない。

　このように、従来の５点法及び８点法によりＥ行列を算出する方法は、６組以上の対応点ペアが与えられた場合に統計的に最適な解を一意に算出することが困難である。

　＜変形例＞
　本実施形態は、上述した例に限定されるものではない。本実施形態は、上述した例に対して、いわゆる当業者が理解し得る多様な変更を適用することが可能である。例えば、本実施形態は、以下の変形例に示す形態によっても実施可能である。

　＜１＞上記の式（７）において自明解を避ける条件は、

に限らない。また、式（７）のその他の式は、すべて任意の定数倍について不変である。そのため、Ｅ行列の１つの成分を１と正規化してもよい。

　＜２＞式（７）の連立方程式を解く際に、式（７）の勾配方向が式（４）のｅ^ＴＭｅを減少させることを条件づけてもよい。例えば、式（７）の勾配方向に対していわゆる直線探索法や信頼領域法を適用して、式（７）がゼロでかつ目的関数を減少させるように、解を探索してもよい。

＜第３実施形態＞
　第３実施形態では、所定の連立方程式を満たすＥ行列のすべての解を算出し、算出したすべての解の中から、複数の対応点ペアについてのエピポーラ誤差が最小となる実数のＥ行列を選択する。

　＜算出装置の構成例＞
　図５は、第３実施形態における算出装置の一例を示すブロック図である。図５において算出装置３０は、取得部１１と、算出部３１とを有している。算出部３１は、局所解算出部３１Ａと、最適解選択部３１Ｂとを有している。

　局所解算出部３１Ａは、取得部１１で取得されたＮ組の対応点ペアを用いて、「所定の連立方程式」を満たすＥ行列のすべての解を算出する。上記の式（７）は、Ｅ行列を未知数とする連立代数方程式である。このため、例えば、グレブナー基底や連続ホモトピー法を用いることによって、上記の式（７）におけるすべての解を計算可能である。

　最適解選択部３１Ｂは、局所解算出部３１Ａで算出されたすべての解の中から、複数の対応点ペアについてのエピポーラ誤差が最小となる実数のＥ行列を１つ選択する。

　＜算出装置の動作例＞
　以上の構成を有する算出装置３０の処理動作の一例について説明する。図６は、第３実施形態における算出装置による処理動作の一例を示すフローチャートである。

　局所解算出部３１Ａは、対応点ペアを用いて、「所定の連立方程式」を満たすＥ行列のすべての解を算出する（ステップＳ２０１）。

　最適解選択部３１Ｂは、ステップＳ２０１で算出されたすべての解の中から、複数の対応点ペアについてのエピポーラ誤差が最小となる実数のＥ行列を選択する（ステップＳ２０２）。

　以上のように第３実施形態によれば、算出装置３０において局所解算出部３１Ａは、「所定の連立方程式」を満たすＥ行列のすべての解を算出する。最適解選択部３１Ｂは、すべての解の中から、複数の対応点ペアについてのエピポーラ誤差が最小となる実数のＥ行列を選択する。

　この算出装置３０の構成により、６組以上の対応点ペアを用いて大域的に最適なＥ行列を一意に算出可能である。その理由は、以下の通りである。例えば、式（７）で表される連立代数方程式の解空間を計算代数幾何に基づいて解析すると、最大で４４０個の実数解を持つことが示される。解析には、例えばMacaulay2などの代数幾何学ソフトウェアを利用することができる。また、グレブナー基底や連続ホモトピー法は、連立方程式の解をすべて列挙できる、ことが知られている。最大で４４０個の実数解のうち、式（４）の目的関数が最小となるＥ行列を１つ選出すれば、そのＥ行列は大域的に最適な唯一の解である。

　＜１＞局所解算出部３１Ａにおいて用いる式（７）の自明解を避ける条件は、

に限らない。例えば、Ｅ行列の１つの成分を１とおいてもよい。その場合、符号の不定性が解消するため、実数解の個数は、最大で２２０個となる。

　＜２＞最適解選択部３１Ｂにおいて大域的最適解を選出する基準となるエピポーラ方程式の誤差は、式（４）の目的関数（代数学的誤差）に限らない。例えば、エピポーラ線と対応点との距離に基づく幾何学的誤差や、幾何学的誤差の１次近似であるサンプソン誤差が用いられてもよい。これらは代数学的誤差よりも厳密な基準として知られており、より高精度なＥ行列を選出できる可能性がある。

＜第４実施形態＞
　第４実施形態は、主に、第１実施形態から第３実施形態のいずれかの算出装置を含むカメラパラメータ推定装置に関する。

　図７は、第４実施形態におけるカメラパラメータ推定装置の一例を示すブロック図である。図７に示すように、カメラパラメータ推定装置５０は、画像供給装置４０と接続されている。カメラパラメータ推定装置５０は、画像供給装置４０と無線で接続されていてもよいし有線で接続されていてもよい。

　画像供給装置４０は、カメラで撮影（形成）された画像を取得して、取得した画像をカメラパラメータ推定装置５０へ供給する。例えば、画像供給装置４０は、ビデオカメラやデジタルカメラを撮影対象となる環境に設置した後に、撮影した画像を画像供給装置４０の不図示の記憶装置（メモリ、ハードディスクなど）に保存し、通信部（不図示）を介してその画像をカメラパラメータ推定装置５０に送信してもよい。

　図７においてカメラパラメータ推定装置５０は、算出装置６０と、分解部５１と、パラメータ算出部５２とを有している。

　算出装置６０は、取得部１１と、算出部６１とを有している。算出装置６０は、第１実施形態から第３実施形態における算出装置１０，２０，３０のいずれに対応していてもよい。すなわち、算出部６１は、算出部１２，２１，３１のいずれに対応していてもよい。

　分解部５１は、算出部６１にて算出されたＥ行列を、並進ベクトルtと回転行列Ｒとに分解する。Ｅ行列の分解方法としては、R. HartleyによるＥ行列の特異値分解を用いる方法や、B. K. P. Hornによる直接解法といった様々な公知の方法を用いてもよい。いずれの方法を用いても、１つのＥ行列は、4組のＲ及びｔへと分解されることが知られている。

　パラメータ算出部５２は、取得部１１からＮ組の対応点ペアを取得する。また、パラメータ算出部５２は、分解部５１から並進ベクトルt及び回転行列Ｒを取得する。そして、パラメータ算出部５２は、並進ベクトルtと回転行列Ｒとを用いて各対応点の３次元座標を三角測量に基づいて計算する。そして、パラメータ算出部５２は、並進ベクトルtと回転行列Ｒと３次元点を出力する。なお、三角測量の方法は広く知られているためここでは説明を省略する。また、４組のR及びtのうち３次元点が２つのカメラの前に復元される１組を選べばよいことも広く知られている。

　＜他の実施形態＞
　図８は、算出装置のハードウェア構成例を示す図である。図８において算出装置１００は、プロセッサ１０１と、メモリ１０２とを有している。プロセッサ１０１は、例えば、マイクロプロセッサ、MPU（Micro Processing Unit）、又はCPU（Central Processing Unit）であってもよい。プロセッサ１０１は、複数のプロセッサを含んでもよい。メモリ１０２は、揮発性メモリ及び不揮発性メモリの組み合わせによって構成される。メモリ１０２は、プロセッサ１０１から離れて配置されたストレージを含んでもよい。この場合、プロセッサ１０１は、図示されていないI/Oインタフェースを介してメモリ１０２にアクセスしてもよい。

　第１実施形態から第３実施形態の算出装置１０，２０，３０は、それぞれ、図８に示したハードウェア構成を有することができる。第１実施形態から第３実施形態の算出装置１０，２０，３０の取得部１１と、算出部１２，２１，３１とは、プロセッサ１０１がメモリ１０２に記憶されたプログラムを読み込んで実行することにより実現されてもよい。プログラムは、様々なタイプの非一時的なコンピュータ可読媒体（non-transitory computer readable medium）を用いて格納され、算出装置１０，２０，３０に供給することができる。非一時的なコンピュータ可読媒体の例は、磁気記録媒体（例えばフレキシブルディスク、磁気テープ、ハードディスクドライブ）、光磁気記録媒体（例えば光磁気ディスク）を含む。さらに、非一時的なコンピュータ可読媒体の例は、ＣＤ－ＲＯＭ（Read Only Memory）、ＣＤ－Ｒ、ＣＤ－Ｒ／Ｗを含む。さらに、非一時的なコンピュータ可読媒体の例は、半導体メモリを含む。半導体メモリは、例えば、マスクＲＯＭ、ＰＲＯＭ（Programmable ROM）、ＥＰＲＯＭ（Erasable PROM）、フラッシュＲＯＭ、ＲＡＭ（Random Access Memory）を含む。また、プログラムは、様々なタイプの一時的なコンピュータ可読媒体（transitory computer readable medium）によって算出装置１０，２０，３０に供給されてもよい。一時的なコンピュータ可読媒体の例は、電気信号、光信号、及び電磁波を含む。一時的なコンピュータ可読媒体は、電線及び光ファイバ等の有線通信路、又は無線通信路を介して、プログラムを算出装置１０，２０，３０に供給できる。

　なお、第４実施形態のカメラパラメータ推定装置５０も、図８に示したハードウェア構成と同様に、プロセッサとメモリとを含むハードウェア構成を有していてもよい。この場合、カメラパラメータ推定装置５０の取得部１１と、算出部６１と、分解部５１と、パラメータ算出部５２とは、プロセッサがメモリに記憶されたプログラムを読み込んで実行することにより実現されてもよい。

　以上、実施の形態を参照して本願発明を説明したが、本願発明は上記によって限定されるものではない。本願発明の構成や詳細には、発明のスコープ内で当業者が理解し得る様々な変更をすることができる。

　１０　算出装置
　１１　取得部
　１２　算出部
　２０　算出装置
　２１　算出部
　２１Ａ　初期値決定部
　２１Ｂ　局所解算出部
　３０　算出装置
　３１　算出部
　３１Ａ　局所解算出部
　３１Ｂ　最適解選択部
　４０　画像供給装置
　５０　カメラパラメータ推定装置
　５１　分解部
　５２　パラメータ算出部
　６０　算出装置
　６１　算出部

Claims

　対応点ペアの２つの対応点間の幾何学的な制約であるエピポーラ制約を表すための基本行列（Ｅ行列）を算出する算出装置であって、
　複数の対応点ペアのそれぞれに含まれる各対応点の画像平面内座標を取得し、各対応点ペアが第１画像及び第２画像にそれぞれ含まれ且つ互いに対応する２つの対応点を含む、取得手段と、
　前記Ｅ行列の特徴を制約条件としたラグランジュ関数Ｌであって前記複数の対応点ペアについてのエピポーラ誤差に対するラグランジュ関数Ｌの１次最適性条件を表し、且つ、ラグランジュ乗数を用いずに前記Ｅ行列の局所解を表す、所定の連立方程式を解くことによって、前記Ｅ行列を算出する算出手段と、
　を具備する算出装置。
　前記算出手段は、行列［Ｍｅ］に対して左から前記Ｅ行列の転置行列を掛けることで得られる行列及び前記行列［Ｍｅ］に対して右から前記Ｅ行列の転置行列を掛けることで得られる行列が対称行列となる第１条件及び前記Ｅ行列の特徴を表す第２条件を含む前記１次最適性条件に基づいて得られる、前記所定の連立方程式を解くことによって、前記Ｅ行列を算出し、
　前記行列［Ｍｅ］は、対応点ペアの第１画像平面内座標及び第２画像平面内座標から算出される係数行列Ｍと前記Ｅ行列を変形したベクトルｅとの積であるベクトルを変形した行列である、
　請求項１記載の算出装置。
　前記算出手段は、前記Ｅ行列の初期値を用いて前記所定の連立方程式を満たすＥ行列の解を算出する局所解算出手段を具備する、
　請求項１又は２に記載の算出装置。
　前記算出手段は、
　前記所定の連立方程式を満たすＥ行列のすべての解を算出する局所解算出手段と、
　前記すべての解の中から、前記複数の対応点ペアについてのエピポーラ誤差が最小となる実数のＥ行列を選択する選択手段と、
　を具備する、
　請求項１又は２に記載の算出装置。
　前記第１画像及び前記第２画像は、内部パラメータが校正済みの１つ以上のカメラによって、同一の物体または同一のシーンが異なる視点から撮影された画像である、
　請求項１から４のいずれか１項に記載の算出装置。
　請求項１から５のいずれか１項に記載の算出装置と、
　前記算出されたＥ行列を、回転行列と並進ベクトルとに分解する分解手段と、
　前記回転行列及び前記並進ベクトルと、各対応点の画像平面内座標とに基づいて、各対応点の３次元座標を算出する３次元座標算出手段と、
　を具備するカメラパラメータ推定装置。
　算出装置によって実行される、対応点ペアの２つの対応点間の幾何学的な制約であるエピポーラ制約を表すための基本行列（Ｅ行列）を算出する算出方法であって、
　複数の対応点ペアのそれぞれに含まれる各対応点の画像平面内座標を取得することを含み、
　各対応点ペアが第１画像及び第２画像にそれぞれ含まれ且つ互いに対応する２つの対応点を含み、
　前記算出方法は、前記Ｅ行列の特徴を制約条件としたラグランジュ関数Ｌであって前記複数の対応点ペアについてのエピポーラ誤差に対するラグランジュ関数Ｌの１次最適性条件を表し、且つ、ラグランジュ乗数を用いずに前記Ｅ行列の局所解を表す、所定の連立方程式を解くことによって、前記Ｅ行列を算出することを含む、
　算出方法。
　対応点ペアの２つの対応点間の幾何学的な制約であるエピポーラ制約を表すための基本行列（Ｅ行列）を算出する処理を、算出装置に実行させるプログラムが格納された非一時的なコンピュータ可読媒体であって、
　前記処理は、複数の対応点ペアのそれぞれに含まれる各対応点の画像平面内座標を取得することを含み、
　各対応点ペアが第１画像及び第２画像にそれぞれ含まれ且つ互いに対応する２つの対応点を含み、
　前記処理は、前記Ｅ行列の特徴を制約条件としたラグランジュ関数Ｌであって前記複数の対応点ペアについてのエピポーラ誤差に対するラグランジュ関数Ｌの１次最適性条件を表し、且つ、ラグランジュ乗数を用いずに前記Ｅ行列の局所解を表す、所定の連立方程式を解くことによって、前記Ｅ行列を算出することを含む、
　非一時的なコンピュータ可読媒体。