特許 6954346 パラメータ推定装置、パラメータ推定方法、及びプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6954346

(24)【登録日】2021年10月4日

(45)【発行日】2021年10月27日

(54)【発明の名称】パラメータ推定装置、パラメータ推定方法、及びプログラム

(51)【国際特許分類】

   G06F 17/17 20060101AFI20211018BHJP

【ＦＩ】

   G06F17/17

【請求項の数】5

【全頁数】18

(21)【出願番号】特願2019-515007(P2019-515007)

(86)(22)【出願日】2017年4月27日

(86)【国際出願番号】JP2017016869

(87)【国際公開番号】WO2018198298

(87)【国際公開日】20181101

【審査請求日】2019年10月21日

(73)【特許権者】

【識別番号】000004237

【氏名又は名称】日本電気株式会社

(74)【代理人】

【識別番号】110002044

【氏名又は名称】特許業務法人ブライタス

(72)【発明者】

【氏名】中野 学

(72)【発明者】

【氏名】柴田 剛志

【審査官】 三坂 敏夫

(56)【参考文献】

【文献】 特表２０１５−５０４６５４（ＪＰ，Ａ）

【文献】 特開２００５−２０１８６９（ＪＰ，Ａ）

【文献】 小柳義夫，ロバスト推定法とデータ解析への応用，オペレーションズ・リサーチ：経営の科学，日本，公益社団法人日本オペレーションズ・リサーチ学会，1978年05月，Vol.23, No.5，第299頁-第304頁，http://www.orsj.or.jp/~archive/menu/03_23.html

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０６Ｆ １７／００−１７／１８

(57)【特許請求の範囲】

【請求項1】

複数のデータ点をアウトライアとインライアとに分離するための閾値を計算して、前記インライアに適合するパラメータを推定するための装置であって、
前記複数のデータ点と前記パラメータの初期値とを取得し、取得した前記初期値に対する前記複数のデータ点それぞれの残差を計算する、残差計算部と、
計算された前記残差に基づいて、前記複数のデータ点それぞれに設定されている重みを更新する、重み更新部と、
前記複数のデータ点に加えて、更新された前記複数のデータ点それぞれの重みも入力として、前記パラメータを推定する、パラメータ推定部と、
計算された前記残差のうちインライアの残差は、期待値と標準偏差の正規分布に従うと仮定し、前記期待値と前記標準偏差と定数に基づいて前記閾値を計算する、閾値決定部と、
推定された前記パラメータと計算された前記閾値とに基づいて、前記パラメータの推定が収束しているかどうかを判定し、収束していないと判定する場合は、推定された前記パラメータを前記初期値として、前記残差計算部、前記重み更新部、前記パラメータ推定部、及び前記閾値決定部それぞれに、再度処理を実行させる、収束判定部と、
を備えている、ことを特徴とするパラメータ推定装置。

【請求項2】

前記収束判定部は、収束していると判定する場合は、推定された前記パラメータを、前記インライアに適合するパラメータとして出力する、
請求項１に記載のパラメータ推定装置。

【請求項3】

前記パラメータ推定部が、反復再重み付け最小二乗法、又はＭ推定を実行して、前記パラメータを推定する、
請求項２に記載のパラメータ推定装置。

【請求項4】

複数のデータ点をアウトライアとインライアとに分離するための閾値を計算して、前記インライアに適合するパラメータを推定するための方法であって、
前記複数のデータ点と前記パラメータの初期値とを取得し、取得した前記初期値に対する前記複数のデータ点それぞれの残差を計算し、
計算された前記残差に基づいて、前記複数のデータ点それぞれに設定されている重みを更新し、
前記複数のデータ点に加えて、更新された前記複数のデータ点それぞれの重みも入力として、前記パラメータを推定し、
計算された前記残差のうちインライアの残差は、期待値と標準偏差の正規分布に従うと仮定し、前記期待値と前記標準偏差と定数に基づいて前記閾値を計算し、
推定された前記パラメータと計算された前記閾値とに基づいて、前記パラメータの推定が収束しているかどうかを判定し、収束していないと判定する場合は、推定された前記パラメータを前記初期値として、前記残差の計算、前記重みの更新、前記パラメータの推定、及び前記閾値の決定それぞれを、再度実行させる、
ことを特徴とするパラメータ推定方法。

【請求項5】

コンピュータによって、複数のデータ点をアウトライアとインライアとに分離するための閾値を計算して、前記インライアに適合するパラメータを推定するためのプログラムであって、
前記コンピュータに、
前記複数のデータ点と前記パラメータの初期値とを取得し、取得した前記初期値に対する前記複数のデータ点それぞれの残差を計算させ、
計算された前記残差に基づいて、前記複数のデータ点それぞれに設定されている重みを更新させ、
前記複数のデータ点に加えて、更新された前記複数のデータ点それぞれの重みも入力として、前記パラメータを推定させ、
計算された前記残差のうちインライアの残差は、期待値と標準偏差の正規分布に従うと仮定し、前記期待値と前記標準偏差と定数に基づいて前記閾値を計算させ、
推定された前記パラメータと計算された前記閾値とに基づいて、前記パラメータの推定が収束しているかどうかを判定し、収束していないと判定する場合は、推定された前記パラメータを前記初期値として、前記残差の計算、前記重みの更新、前記パラメータの推定、及び前記閾値の決定それぞれを、再度実行させる、
ことを特徴とするプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、外れ値を含む入力データ点に適合する幾何パラメータを推定するための、パラメータ推定装置、及びパラメータ推定方法に関し、更には、これらを実現するためのプログラムに関する。

【背景技術】

【0002】

複数のデータ点に適合する幾何パラメータを推定することは、回帰分析、又はパラメータフィッティングなどと呼ばれ、統計学、信号処理、画像処理、コンピュータビジョンといった幅広い分野における重要な要素技術となっている。推定するパラメータとしては、例えば、データ点が2次元であれば直線、放物線、ロジスティック曲線等が挙げられ、データ点が3次元であれば平面、曲面、球等が挙げられ、更に、データ点が多次元であれば超平面等が挙げられる。以下では、パラメータの種類は特に限定せず、一般的な場合について説明する。

【0003】

入力となるデータ点は一般に観測ノイズを含むため、求めるべきパラメータ上の値とは完全には適合しない。これは、パラメータとデータ点との間に誤差（または「距離」とも呼ばれる）が生じているからである。そして、観測ノイズが期待値0の正規分布に従うと予想されるとき、全データ点の誤差を最小にするパラメータを推定する方法としては、最小二乗法が広く用いられている。最小二乗法が用いられる理由は、データ点の数が多ければ、観測ノイズは正規分布に近似できるという中心極限定理が成立すると考えられるからである。なお、分布が既知である限定された状況においては、データ点の事前に近似できる分布は正規分布に限らないが、以下では正規分布として説明する。

【0004】

実際のパラメータ推定においては、データ点の中で、他のデータ点から大きく外れた値を持つアウトライア、（outlier、外れ値とも呼ばれる）が観測されることがある。一方で、外れ値ではないデータ点、つまり、それが含む観測ノイズが正規分布に従っており、パラメータ上の値に正しく適合する点を、インライア（inlier）と呼ぶ。アウトライアが生じる原因は、例えば、観測ミスや記録ミス、または正規分布に当てはまらない特異な個体などが原因である。最小二乗法は、全データ点が正規分布という前提のため、このようなアウトライアが1点でもあると、その大きな誤差を最小化しようと働く。また、この結果、その他のインライアに対して当てはまりの精度が低いパラメータが推定されてしまう。

【0005】

このようなアウトライアを除去してパラメータを高精度に推定する方法として、非特許文献１に開示されているLMeds (Least Median Squares)、非特許文献２に開示されているRANSAC (Random Sample Consensus)、非特許文献３に開示されているM推定（M-estimator）といった推定方法が広く知られている。

【0006】

非特許文献１に開示されているLMedsは、データ点からパラメータ推定に必要な最小点数をランダムに抽出して、パラメータ推定を行い、全点の誤差と誤差の中央値とを計算する。そして、LMedsは、ランダムサンプリングを一定数繰り返した後、誤差の中央値が最小となるパラメータを正しい推定とみなす。LMedsによれば、アウトライアが最大で50%含まれていても対応可能である。

【0007】

また、非特許文献２に開示されているRANSACも、LMedsと同様に、データ点からパラメータ推定に必要な最小点数をランダムに抽出する。但し、RANSACは、事前に決定した閾値以下の誤差であるインライアを数え上げる。そして、RANSACは、ランダムサンプリングとインライアの数え上げとを一定数繰り返した後、インライアの数が最大となるパラメータを正しい推定とみなす。RANSACは、アウトライアが50%以上であっても対応でき、計算量がLMedsよりも少ないのが利点である。

【0008】

M推定は、最小二乗法における二乗誤差の代わりに、データ点の誤差が大きいほど寄与率が小さくなるような重み付け関数を用いて、誤差の計算と重みの更新とを、誤差が収束するまで反復する。なお、非特許文献4によると、M推定はIRLS （Iterative Reweighted Least Squares：反復再重み付け最小二乗法）と数学的に等価なことが示されている。重み付けの大きさを制御するための閾値としては、M推定でも、LMeds及びRANSACのように固定値が用いられる。つまり、誤差が閾値より大きければ重みを小さくし、閾値より小さければ重みを大きくする。IRLSでは、閾値を十分大きな値で初期化し、反復のたびに0より大きく1より小さい減衰定数を掛けて徐々に閾値を減少させる。M推定は、IRLSにおいて減衰定数を1に固定した場合に等しいため、以下において、M推定の代わりにIRLSについて説明する。

【先行技術文献】

【非特許文献】

【0009】

【非特許文献1】P. J Rousseeuw, "Least Median of Squares Regression". Journal of the American Statistical Association. 79: 871-880, 1984.

【非特許文献2】M. A. Fischler and R. C. Bolles, "Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography," Comm. of the ACM, vol.24, no.6, pp.381-395, 1981.

【非特許文献3】Z. Zhang, “Parameter estimation techniques: A tutorial with application to conic fitting,” International journal of Image and Vision Computing, vol.15, no.1, pp.59-76, 1997.

【非特許文献4】M. Black and A. Rangarajan, “On the unification of line processes, outlier rejection, and robust statistics with applications in early vision,” International journal of computer vision, vol.19, no.1, pp.57-91, 1996.

【発明の概要】

【発明が解決しようとする課題】

【0010】

しかしながら、LMeds及びRANSACといったランダムサンプリングに基づく手法には、以下のような問題がある。以下に具体的に説明する。

【0011】

まず、必要な繰り返し数がアウトライアの割合に依存して指数関数的に増大する。そのため、新規のデータに対して事前に計算量を見積もることが難しい。特に、アウトライア率が50%を超えるような場合、LMedsは原理的に対応できず、RANSACは対応できるが常に高速性が要求される応用には適さない。

【0012】

次に、ランダムサンプリングは原理的に同一の結果を得られず、同じデータ点であっても複数回実行した場合は結果が異なるため、安定性が低い。そして、通常はアウトライアの誤差分布が不明なため、閾値の適切な決定が困難である。閾値が大きすぎるとアウトライアをインライアと誤判定してしまう。小さすぎるとインライアをアウトライアと誤判定してしまう上に、反復回数がさらに増加する一因となる。

【0013】

一方、IRLSの反復回数は、閾値の初期値と、減衰定数と、事前に定義した最小の閾値と、により決定されるため、アウトライアの割合に関係なく最大の計算量を見積もることは可能である。また、ランダムサンプリングを行わないため、、IRLSは、同一のデータ点に対しては同一の結果を返す。このように、IRLSを用いれば、上述したLMedsとRANSACとの問題の一部は解消される。

【0014】

但し、IRLSを用いた場合であっても、アウトライアの誤差分布は不明であるため、閾値と減衰定数とを適切に決定することは困難である。更に、IRLSにおいては、反復毎に閾値が連続的に変更されていくため、LMedsとRANSACとの場合よりも、閾値の影響はより大きなものとなる。

【0015】

本発明の目的の一例は、上記問題を解消し、複数のデータ点に適合する幾何パラメータの推定において、アウトライアとインライアとを分離する際の閾値を適切に決定し得る、パラメータ推定装置、パラメータ推定方法、およびプログラムを提供することである。

【課題を解決するための手段】

【0016】

上記目的を達成するため、本発明の一側面におけるパラメータ推定装置は、
複数のデータ点をアウトライアとインライアとに分離するための閾値を計算して、前記インライアに適合するパラメータを推定するための装置であって、
前記複数のデータ点を入力として、前記パラメータを推定する、パラメータ推定部と、
前記複数のデータ点それぞれの残差の統計情報に基づいて、前記閾値を計算する、閾値決定部と、
推定された前記パラメータと計算された前記閾値とに基づいて、前記パラメータの推定が収束しているかどうかを判定し、収束していないと判定する場合は、前記パラメータ推定部、及び前記閾値決定部それぞれに、再度処理を実行させる、収束判定部と、
を備えている、ことを特徴とする。

【0017】

また、上記目的を達成するため、本発明の一側面におけるパラメータ推定方法は、
複数のデータ点をアウトライアとインライアとに分離するための閾値を計算して、前記インライアに適合するパラメータを推定するための方法であって、
（ａ）前記複数のデータ点を入力として、前記パラメータを推定する、ステップと、
（ｂ）前記複数のデータ点それぞれの残差の統計情報に基づいて、前記閾値を計算する、ステップと、
（ｃ）推定された前記パラメータと計算された前記閾値とに基づいて、前記パラメータの推定が収束しているかどうかを判定する、ステップと、を有し、
前記（ｃ）のステップで、収束していないと判定する場合は、前記（ａ）のステップ、及び前記（ｂ）のステップそれぞれが再度実行される、
ことを特徴とする。

【0018】

更に、上記目的を達成するため、本発明の一側面におけるプログラムは、
コンピュータによって、複数のデータ点をアウトライアとインライアとに分離するための閾値を計算して、前記インライアに適合するパラメータを推定するためのプログラムであって、
前記コンピュータに、
（ａ）前記複数のデータ点を入力として、前記パラメータを推定する、ステップと、
（ｂ）前記複数のデータ点それぞれの残差の統計情報に基づいて、前記閾値を計算する、ステップと、
（ｃ）推定された前記パラメータと計算された前記閾値とに基づいて、前記パラメータの推定が収束しているかどうかを判定する、ステップと、
を実行させ、
前記（ｃ）のステップで、収束していないと判定する場合は、前記（ａ）のステップ、及び前記（ｂ）のステップそれぞれを再度実行する、
ことを特徴とする。

【発明の効果】

【0019】

以上のように本発明によれば、複数のデータ点に適合する幾何パラメータの推定において、アウトライアとインライアとを分離する際の閾値を適切に決定することができる。

【図面の簡単な説明】

【0020】

【図1】図１は、本発明の実施の形態のけるパラメータ推定装置の概略構成を示すブロック図である。

【図2】図２は、本発明の実施の形態におけるパラメータ推定装置の具体的構成を示すブロック図である。

【図3】図３は、本発明の実施の形態におけるパラメータ推定装置の動作を示すフロー図である。

【図4】図４は、本発明の実施の形態におけるインライアとアウトライアとの残差分布の遷移を示す図である。

【図5】図５は、本発明の実施の形態におけるパラメータ推定装置を実現するコンピュータの一例を示すブロック図である。

【発明を実施するための形態】

【0021】

（実施の形態）
以下、本発明の実施の形態における、パラメータ推定装置、パラメータ推定方法、プログラムについて、図１〜図５を参照しながら説明する。

【0022】

［装置構成］
最初に、図１を用いて、本実施の形態におけるパラメータ推定装置１０の構成について説明する。図１は、本発明の実施の形態のけるパラメータ推定装置の概略構成を示すブロック図である。

【0023】

図１に示す、本実施の形態におけるパラメータ推定装置１０は、複数のデータ点をアウトライアとインライアとに分離するための閾値を計算して、インライアに適合するパラメータを推定する装置である。図１に示すように、パラメータ推定装置１０は、パラメータ推定部１３と、閾値決定部１４と、収束判定部１５とを備えている。

【0024】

パラメータ推定部１３は、複数のデータ点を入力として、パラメータを推定する。閾値決定部１４は、複数のデータ点それぞれの残差の統計情報に基づいて、閾値を計算する。収束判定部１５は、推定されたパラメータと計算された閾値とに基づいて、パラメータの推定が収束しているかどうかを判定する。また、収束判定部１５は、収束していないと判定する場合は、パラメータ推定部１３、及び閾値決定部１４それぞれに、再度処理を実行させる。

【0025】

このように、本実施の形態では、複数のデータ点をアウトライアとインライアとに分離するための閾値は、パラメータを反復して推定する際において、残差の統計情報に基づいて計算されている。このため、本実施形態によれば、複数のデータ点に適合する幾何パラメータの推定において、アウトライアとインライアとを分離する際の閾値、具体的には、IRLSにおける閾値を適切に自動的に決定することができる。

【0026】

続いて、図２を用いて、本実施の形態におけるパラメータ推定装置１０の構成をより具体的に説明する。図２は、本発明の実施の形態におけるパラメータ推定装置の具体的構成を示すブロック図である。図２に示すように、本実施の形態では、パラメータ推定装置１０は、パラメータ推定部１３、閾値決定部１４、及び収束判定部１５に加えて、残差計算部１１と、重み更新部１２とを更に備えている。

【0027】

残差計算部１１は、まず、入力として、複数のデータ点とパラメータの初期値とを取得する。また、残差計算部１１は、取得した初期値に対する複数のデータ点それぞれの残差を計算する。

【0028】

重み更新部１２は、計算された残差に基づいて、複数のデータ点それぞれに設定されている重みを更新する。また、本実施の形態では、重み更新部１２は、事前に決定された重み関数を用いて、各データ点の重みを更新する。

【0029】

パラメータ推定部１３は、本実施の形態では、複数のデータ点に加えて、更新された複数のデータ点それぞれの重みも入力として、重み付のデータ点に適合するパラメータを推定する。

【0030】

閾値決定部１４は、本実施の形態では、残差計算部１１が計算した残差のうちインライアの残差を特定し、特定したインライアの残差の統計情報に基づいて、閾値を計算する。具体的には、閾値決定部１４は、まず、閾値の初期値を取得する。次いで、閾値決定部１４は、取得した初期値を、統計情報を用いて更新して、新たな閾値を計算する。その後、再度の閾値の計算が必要となると、閾値決定部１４は、前回更新された閾値を更に更新する。

【0031】

また、本実施の形態において、統計情報としては、インライアの残差の平均値、標準偏差、中央値、歪度、尖度、エントロピー等が挙げられる。統計情報は、インライアの残差が持つと予想される確率分布に応じて、異なる値となる。

【0032】

収束判定部１５は、本実施の形態では、収束していると判定する場合は、推定されたパラメータを、インライアに適合するパラメータとして外部に出力する。また、収束判定部１５は、収束していないと判定する場合は、推定されたパラメータを初期値として、残差計算部１１、重み更新部１２、パラメータ推定部１３、及び閾値決定部１４それぞれに、再度処理を実行させる。つまり、本実施の形態では、パラメータの推定が収束するまで、パラメータの初期値を更新しながら、反復してパラメータの推定が行なわれる。なお、収束判定部１５における判定処理の具体例については後述する。

【0033】

このように、本実施の形態において、パラメータ推定装置１０は、データ点をインライアとアウトライアに分離する閾値をパラメータに対するデータ点の残差に基づき決定し、更に、決定を行なう処理を、残差が収束するまで反復する。そして、収束すると、パラメータ推定装置１０は、インライアに適合するパラメータを出力する。

【0034】

［装置構成］
次に、本実施の形態におけるパラメータ推定装置１０の動作について図３を用いて説明する。図３は、本発明の実施の形態におけるパラメータ推定装置の動作を示すフロー図である。以下の説明においては、適宜図１及び図２を参酌する。また、本実施の形態では、パラメータ推定装置１０を動作させることによって、パラメータ推定方法が実施される。よって、本実施の形態におけるパラメータ推定方法の説明は、以下のパラメータ推定装置１０の動作説明に代える。

【0035】

図３に示すように、最初に、パラメータ推定装置１０に、複数のデータ点と、パラメータの初期値と、閾値の初期とが入力されると、残差計算部１１は、データ点毎に、パラメータの初期値に対する残差を計算する（ステップＳ１）。

【0036】

ステップＳ１において、パラメータの初期値としては、例えば、乱数が与えられても良いし、事前に決定されているデフォルト値（例えば、直線であれば、「傾き１」、「切片０」）が与えられても良い。また、データ点とパラメータとが時系列により連続的に変化するときは、前の時刻において計算されたパラメータが初期値として入力されてもよい。また、画像特徴点のマッチングスコアように、データ点の信頼度が事前に得られている場合は、信頼度による重み付けを用いた最小二乗法によって初期値が与えられても良い。

【0037】

また、閾値の初期値としては、例えば、プログラムにおいて取りうる最大の実数値が与えられても良いし、残差の上限が実験的または経験的に知られている場合は、それらが与えられていても良い。また、データ点とパラメータとが時系列により連続的に変化するときは、前の時刻における残差の最大値を定数倍して得られた値が与えられていても良い。

【0038】

更に、入力されるデータ点の数は、推定すべきパラメータを表現するために必要最小限の数以上であれば良い。例えば、パラメータが直線であれば２点以上、平面であれば３点以上である。また、もし点の数がパラメータの推定に足りない場合は、パラメータ推定装置１０は、実行不可能であることを示すフラグ、または実行不可能であることを示すパラメータ（例えば、全ての値がゼロまたは無限大）を出力して、処理を終了しても良い。

【0039】

また、本実施の形態において、残差は、パラメータに対する各データ点の適合度合いを表す指標である。残差としては、例えば、データ点からパラメータが表す線又は面に下ろした垂線がなす幾何学的距離、データ点のノイズ分散を考慮したマハラノビス距離等が挙げられる。

【0040】

次いで、重み更新部１２は、ステップＳ１で計算された残差を入力として、事前に決定されている重み関数に従って、各データ点の重みを更新する（ステップＳ２）。ステップＳ２で用いられる重み関数としては、残差に反比例した値を返却する関数、具体的には、残差が大きいほど小さい値を返し、残差が小さいほど大きい値を返す関数が挙げられる。

【0041】

次いで、パラメータ推定部１３は、ステップＳ２で更新された各データ点の重みと、データ点とを入力として、重み付きデータ点に適合するパラメータを推定する（ステップＳ３）。

【0042】

次いで、閾値決定部１４は、ステップＳ１で計算された残差のうち、インライアの残差を入力として、その統計情報を利用して、新たな閾値を決定する（ステップＳ４）。

【0043】

次に、収束判定部１５は、ステップＳ３で推定されたパラメータと、ステップＳ４で計算された閾値とに基づいて、パラメータの推定が収束しているかどうかを判定する（ステップＳ５）。

【0044】

ステップＳ５の判定の結果、収束していない場合は、収束判定部１５は、ステップＳ３で推定されたパラメータをパラメータの初期値とし、ステップＳ４で計算された閾値を閾値の初期値とする。その後、残差計算部１１が、再度、新たな初期値を用いて、再度ステップＳ１を実行する。これにより、ステップＳ２〜Ｓ４も再度実行される。

【0045】

一方、ステップＳ５の判定の結果、収束している場合は、収束判定部１５は、ステップＳ３で推定されたパラメータを出力する（ステップＳ６）。これにより、パラメータ推定装置１０における処理は終了する。

【0046】

また、本実施の形態における「パラメータの推定の収束」とは、一般の最適化計算又は数値計算において用いられる「収束」と同様の意味であり、これ以上反復計算を繰り返してもパラメータ及び閾値の値が変化しない、と判断される状態にあることをいう。

【0047】

従って、収束判定部１５は、例えば、一つ前のステップＳ１〜Ｓ４までの処理と最新のステップＳ１〜Ｓ４までの処理とにおける、パラメータの差分量及び閾値の差分量のうちの一方又は両方を計算し、差分量が一定値以下であるときに収束したと判定することができる。また、収束判定部１５は、例えば、最新のステップＳ４で決定された閾値が一定値以下であるときに、収束したと判定することもできる。

【0048】

［具体例］
続いて、パラメータ推定装置１０における処理の具体例について以下に説明する。なお、以下の具体例においては、推定するパラメータをベクトルθ、i番目のデータ点をベクトルx、パラメータθに対するi番目のデータ点の残差を計算する関数をf(x_i;θ)とする。また、i番目のデータ点の重みをw_i、閾値をλとし、データ点は全部でN点あるとする。なお、重みは各データ点の間で相対的なものであるため、本具体例では、0から1の間の実数（0≦ w_i ≦1）とする。

【0049】

まず、IRLSの定義を行う。IRLSは、以下の数１で表される最適化を、閾値λを初期値から徐々に小さくしながら繰り返す操作のことである。

【0050】

【数1】

【0051】

ここで、関数g(w_i)は、penalty functionと呼ばれ、w_iが1に近づくと0、1に近づくと0となる任意の関数である。本具体例では、以下の数２が成立するとして説明する。

【0052】

【数2】

【0053】

ある反復（ステップＳ１〜Ｓ４）において、各重みw_iは、上記数１をw_iについて微分してゼロとした以下の数３で表される方程式を解くことで得られる。

【0054】

【数3】

【0055】

これを上記数２の場合について解くと、以下の数４が得られる。

【0056】

【数4】

【0057】

上記数４をλの関数と見れば、0から1の値をとることが確認できる。また、上記数４はGeman-McClureによるM推定の重み付け関数と一致することが、上述した非特許文献４に示されている。

【0058】

ここで、上述した具体例を、図３に示したステップＳ１〜Ｓ５に沿って更に説明する。まず、ステップＳ１において、残差計算部１１は、パラメータθに対する各データ点x_iの残差f(x_i;θ)を計算する。

【0059】

次に、ステップＳ２において、重み更新部１２は、閾値λと、上記数４に従って、各データ点x_iの重みw_iを更新する。

【0060】

次に、ステップＳ３において、パラメータ推定部１３は、ステップＳ２で更新された重みw_iを上記数１に代入して、最適化問題を最小化するパラメータθを推定する。最適化問題の解は、例えば、上記数１がｍf(x_i;θ)の二乗和であることを利用して、ガウスニュートン法、レーベンバーグマーカード法などによって求めることができる。

【0061】

次に、ステップＳ４において、閾値決定部１４は、ステップＳ１で計算された残差f(x_i;θ)と、ステップＳ２で更新された重みw_iとを用いて、新たな閾値λを決定する。

【0062】

具体的には、閾値決定部１４は、まず、重みw_iから各データ点がインライアかどうかを判定する。重みは0から1の間の値を取るから、閾値決定部１４は、例えば重みw_iが0.5以上のデータ点をインライアとして判定する。なお、重みw_iが１に近い値となるときにインライアとして判定するとすれば、より厳しい判定結果となる。

【0063】

次に、インライアの残差だけを用いて、閾値決定部１４は、新たなしきい値λを決定する。ここで、インライアの残差が、期待値μ、標準偏差σの正規分布に従うと仮定すると、インライアの残差は、ある定数βを用いてμ＋βσの範囲内に分布すると考えられる。例えば、β=1、β=2とすると、残差f(x_i;θ)は、それぞれ約84％、約97%となり、これらはμ＋βσ以下となる。つまり、λ=μ＋βσと更新すれば、その反復時点の任意の割合のインライアを抽出できる。

【0064】

次に、ステップＳ５において、収束判定部１５は、ステップＳ３で推定されたパラメータθとステップＳ４で決定された閾値λとを入力として、パラメータ推定が収束したかどうかを判定する。

【0065】

上述したように、例えば、収束判定部１５、前回のステップＳ１〜Ｓ４の実行によって得られたパラメータθ又は閾値λと、前回のステップＳ１〜Ｓ４の実行によって得られたパラメータθ又は閾値λとの差分を比較し、差分が一定値以下であれば収束したと判定する。また、収束判定部１５は、最新のステップＳ４で決定された閾値λが事前に決めていた最小値よりも小さくなった場合に、収束したと判定しても良い。

【0066】

判定の結果、収束していない場合は、収束判定部１５は、ステップＳ３で推定されたパラメータをパラメータの初期値とし、ステップＳ４で計算された閾値を閾値の初期値とする。その後、再度ステップＳ１以降が実行される。一方、収束している場合は、収束判定部１５は、ステップＳ３で推定されたパラメータを出力する

【0067】

［実施の形態における効果］
以上のように、本実施の形態によれば、IRLSにおける閾値を自動的に且つ適切に決定でき、アウトライアとインライアとを高精度に分離できる。また、この結果、インライアに適合するパラメータを高精度に推定可能である。

【0068】

ここで、上述したステップＳ１〜Ｓ４を反復して実行することにより、インライアとアウトライアとの分離ができる原理を、図４を用いて説明する。図４は、本発明の実施の形態におけるインライアとアウトライアとの残差分布の遷移を示す図である。図４は、上段から下段へと初期状態から収束状態までを示している。

【0069】

図４に示すように、上段の初期状態、つまり、インライアとアウトライアの区別がほとんどつかない状態では、２つの残差分布が混じっている。中心極限定理により、これは混合正規分布と仮定できる。ここで、ステップＳ４においてλ=μ＋βσとすることは、μ＋βσより大きい残差f(x_i;θ)を持つデータ点の重みw_iが急激に小さくなることを意味する。そして、データ点により大きな重みが与えられることは、パラメータにおいてインライアの占める割合が上昇することと等価である。インライア率が高まった状態で、上記の数１を解くと、よりインライアに適合するパラメータθが求まり、図４の中段の状態となる。そして、ステップＳ１〜Ｓ４が反復されると、最終的に図４の下段に示すように、インライアとアウトライアとが分離された状態に収束する。

【0070】

［変形例］
ここで、本実施の形態における変形例１〜４について説明する。本発明は、上述した実施の形態のみに限定されることはない。本発明は、上述した実施の形態に対して、いわゆる当業者が理解し得る多様な変更を適用することが可能である。例えば、本発明は、以下の変形例に示す形態によっても実施可能である。

【0071】

（１）変形例１
どのpenalty functionを用いるかは、上述した実施の形態に限定されない。例えば、非特許文献４にはpenalty functionの様々な具体例が詳しく記載されている。また、非特許文献４にはIRLSとM推定との間の変換方法が記載されているため、非特許文献3に記載されている様々なM推定の関数が利用できる。

【0072】

（２）変形例２
閾値決定方法は、上述した実施の形態に限定されない。例えば、前回の反復（ステップＳ１〜Ｓ４の実施）で得られた閾値を保存しておき、この閾値と今回得られた閾値との差分が一定値以下となるようにして、急激な減衰が抑制されても良い。また、例えば、従来の減衰定数を用いる方法も同時に計算しておき、より小さな値又はより大きな値が、新たな閾値として用いられても良い。また、例えば、インライアの残差分布がコーシー分布である場合、期待値と標準偏差とは存在しないため、代わりに最頻値と半値半幅とを与えるスケール定数が用いられても良い。

【0073】

（３）変形例３
収束判定方法は、上述した実施の形態に限定されない。例えば、収束の判定基準として、反復の前後における上記数１で表される目的関数の差分が用いられても良いし、データ点の重みw_iの差分が用いられても良い。このとき、収束判定部１５は、上記数１で表される目的関数または重みw_iを入力として受け付ける。

【0074】

（４）変形例４
パラメータ推定装置１０に入力されるデータは、上述の実施の形態に示された例に限定されない。例えば、パラメータの初期値の代わりに、重みの初期値が入力されても良い。この場合、パラメータ推定装置１０においては、パラメータ推定部１３が最初に処理を実行することになり、次いで残差計算部１１、重み更新部１２、閾値決定部１４、収束判定部１５の順で処理が実行される。なお、各部による処理の実行順序は、適宜変更可能である。

【0075】

［プログラム］
本実施の形態におけるプログラムは、コンピュータに、図３に示すステップＳ１〜Ｓ６を実行させるプログラムであれば良い。このプログラムをコンピュータにインストールし、実行することによって、本実施の形態におけるパラメータ装置とパラメータ推定方法とを実現することができる。この場合、コンピュータのＣＰＵ（Central Processing Unit）は、残差計算部１１、重み更新部１２、パラメータ推定部１３、閾値決定部１４、及び収束判定部１５として機能し、処理を行なう。

【0076】

また、本実施の形態におけるプログラムは、複数のコンピュータによって構築されたコンピュータシステムによって実行されても良い。この場合は、例えば、各コンピュータが、それぞれ、残差計算部１１、重み更新部１２、パラメータ推定部１３、閾値決定部１４、及び収束判定部１５のいずれかとして機能しても良い。

【0077】

ここで、本実施の形態におけるプログラムを実行することによって、パラメータ推定装置１０を実現するコンピュータについて図５を用いて説明する。図５は、本発明の実施の形態におけるパラメータ推定装置を実現するコンピュータの一例を示すブロック図である。

【0078】

図５に示すように、コンピュータ１１０は、ＣＰＵ１１１と、メインメモリ１１２と、記憶装置１１３と、入力インターフェイス１１４と、表示コントローラ１１５と、データリーダ／ライタ１１６と、通信インターフェイス１１７とを備える。これらの各部は、バス１２１を介して、互いにデータ通信可能に接続される。

【0079】

ＣＰＵ１１１は、記憶装置１１３に格納された、本実施の形態におけるプログラム（コード）をメインメモリ１１２に展開し、これらを所定順序で実行することにより、各種の演算を実施する。メインメモリ１１２は、典型的には、ＤＲＡＭ（Dynamic Random Access Memory）等の揮発性の記憶装置である。また、本実施の形態におけるプログラムは、コンピュータ読み取り可能な記録媒体１２０に格納された状態で提供される。なお、本実施の形態におけるプログラムは、通信インターフェイス１１７を介して接続されたインターネット上で流通するものであっても良い。

【0080】

また、記憶装置１１３の具体例としては、ハードディスクドライブの他、フラッシュメモリ等の半導体記憶装置が挙げられる。入力インターフェイス１１４は、ＣＰＵ１１１と、キーボード及びマウスといった入力機器１１８との間のデータ伝送を仲介する。表示コントローラ１１５は、ディスプレイ装置１１９と接続され、ディスプレイ装置１１９での表示を制御する。

【0081】

データリーダ／ライタ１１６は、ＣＰＵ１１１と記録媒体１２０との間のデータ伝送を仲介し、記録媒体１２０からのプログラムの読み出し、及びコンピュータ１１０における処理結果の記録媒体１２０への書き込みを実行する。通信インターフェイス１１７は、ＣＰＵ１１１と、他のコンピュータとの間のデータ伝送を仲介する。

【0082】

また、記録媒体１２０の具体例としては、ＣＦ（Compact Flash（登録商標））及びＳＤ（Secure Digital）等の汎用的な半導体記憶デバイス、フレキシブルディスク（Flexible Disk）等の磁気記録媒体、又はＣＤ−ＲＯＭ（Compact Disk Read Only Memory）などの光学記録媒体が挙げられる。

【0083】

なお、本実施の形態におけるパラメータ推定装置１０は、プログラムがインストールされたコンピュータではなく、各部に対応したハードウェアを用いることによっても実現可能である。更に、パラメータ推定装置１０は、一部がプログラムで実現され、残りの部分がハードウェアで実現されていてもよい。

【0084】

上述した実施の形態の一部又は全部は、以下に記載する（付記１）〜（付記１５）によって表現することができるが、以下の記載に限定されるものではない。

【0085】

（付記１）
複数のデータ点をアウトライアとインライアとに分離するための閾値を計算して、前記インライアに適合するパラメータを推定するための装置であって、
前記複数のデータ点を入力として、前記パラメータを推定する、パラメータ推定部と、
前記複数のデータ点それぞれの残差の統計情報に基づいて、前記閾値を計算する、閾値決定部と、
推定された前記パラメータと計算された前記閾値とに基づいて、前記パラメータの推定が収束しているかどうかを判定し、収束していないと判定する場合は、前記パラメータ推定部、及び前記閾値決定部それぞれに、再度処理を実行させる、収束判定部と、
を備えている、ことを特徴とするパラメータ推定装置。

【0086】

（付記２）
前記複数のデータ点と前記パラメータの初期値とを取得し、取得した前記初期値に対する前記複数のデータ点それぞれの残差を計算する、残差計算部と、
計算された前記残差に基づいて、前記複数のデータ点それぞれに設定されている重みを更新する、重み更新部と、を更に備え、
前記パラメータ推定部は、前記複数のデータ点に加えて、更新された前記複数のデータ点それぞれの重みも入力として、前記パラメータを推定し、
前記閾値決定部は、計算された前記残差のうち前記インライアの残差を特定し、特定した前記インライアの残差の前記統計情報に基づいて、前記閾値を計算し、
前記収束判定部は、
収束していると判定する場合は、推定された前記パラメータを、前記インライアに適合するパラメータとして出力し、
収束していないと判定する場合は、推定された前記パラメータを前記初期値として、前記残差計算部、前記重み更新部、前記パラメータ推定部、及び前記閾値決定部それぞれに、再度処理を実行させる、
付記１に記載のパラメータ推定装置。

【0087】

（付記３）
前記パラメータ推定部が、反復再重み付け最小二乗法、又はM推定を実行して、前記パラメータを推定する、
付記２に記載のパラメータ推定装置。

【0088】

（付記４）
前記閾値決定部は、前記インライアの残差の事前確率分布から前記統計情報を求め、求めた前記統計情報に基づいて、前記閾値を計算する、
付記２または３に記載のパラメータ推定装置。

【0089】

（付記５）
前記事前確率分布は、混合正規分布で表され、
前記閾値決定部は、前記統計情報として、前記インライアの残差の平均値と標準偏差とを求める、
付記４に記載のパラメータ推定装置。

【0090】

（付記６）
複数のデータ点をアウトライアとインライアとに分離するための閾値を計算して、前記インライアに適合するパラメータを推定するための方法であって、
（ａ）前記複数のデータ点を入力として、前記パラメータを推定する、ステップと、
（ｂ）前記複数のデータ点それぞれの残差の統計情報に基づいて、前記閾値を計算する、ステップと、
（ｃ）推定された前記パラメータと計算された前記閾値とに基づいて、前記パラメータの推定が収束しているかどうかを判定する、ステップと、を有し、
前記（ｃ）のステップで、収束していないと判定する場合は、前記（ａ）のステップ、及び前記（ｂ）のステップそれぞれが再度実行される、
ことを特徴とするパラメータ推定方法。

【0091】

（付記７）
（ｄ）前記複数のデータ点と前記パラメータの初期値とを取得し、取得した前記初期値に対する前記複数のデータ点それぞれの残差を計算する、ステップと、
（ｅ）計算された前記残差に基づいて、前記複数のデータ点それぞれに設定されている重みを更新する、ステップと、を更に有し、
前記（ａ）のステップで、前記複数のデータ点に加えて、前記（ｅ）のステップで更新された前記複数のデータ点それぞれの重みも入力として、前記パラメータを推定し、
前記（ｂ）のステップで、前記（ｄ）のステップで計算された前記残差のうち前記インライアの残差を特定し、特定した前記インライアの残差の前記統計情報に基づいて、前記閾値を計算し、
前記（ｃ）のステップで、収束していると判定する場合は、前記（ｃ）のステップにおいて、推定された前記パラメータを、前記インライアに適合するパラメータとして出力し、
前記（ｃ）のステップで、収束していないと判定する場合は、前記（ａ）のステップで推定された前記パラメータを前記初期値として、前記（ｄ）のステップ、前記（ｅ）のステップ、前記（ａ）のステップ、及び前記（ｂ）のステップそれぞれが再度実行される、
付記６に記載のパラメータ推定方法。

【0092】

（付記８）
前記（ａ）のステップにおいて、反復再重み付け最小二乗法、又はM推定を実行して、前記パラメータを推定する、
付記７に記載のパラメータ推定方法。

【0093】

（付記９）
前記（ｂ）のステップにおいて、前記インライアの残差の事前確率分布から前記統計情報を求め、求めた前記統計情報に基づいて、前記閾値を計算する、
付記７または８に記載のパラメータ推定方法。

【0094】

（付記１０）
前記事前確率分布は、混合正規分布で表され、
前記（ｂ）のステップにおいて、前記統計情報として、前記インライアの残差の平均値と標準偏差とを求める、
付記９に記載のパラメータ推定方法。

【0095】

（付記１１）
コンピュータによって、複数のデータ点をアウトライアとインライアとに分離するための閾値を計算して、前記インライアに適合するパラメータを推定するためのプログラムであって、
前記コンピュータに、
（ａ）前記複数のデータ点を入力として、前記パラメータを推定する、ステップと、
（ｂ）前記複数のデータ点それぞれの残差の統計情報に基づいて、前記閾値を計算する、ステップと、
（ｃ）推定された前記パラメータと計算された前記閾値とに基づいて、前記パラメータの推定が収束しているかどうかを判定する、ステップと、
を実行させ、
前記（ｃ）のステップで、収束していないと判定する場合は、前記（ａ）のステップ、及び前記（ｂ）のステップそれぞれを再度実行する、
ことを特徴とするプログラム。

【0096】

（付記１２）
前記コンピュータに、
（ｄ）前記複数のデータ点と前記パラメータの初期値とを取得し、取得した前記初期値に対する前記複数のデータ点それぞれの残差を計算する、ステップと、
（ｅ）計算された前記残差に基づいて、前記複数のデータ点それぞれに設定されている重みを更新する、ステップと、
を更に実行させ、
前記（ａ）のステップで、前記複数のデータ点に加えて、前記（ｅ）のステップで更新された前記複数のデータ点それぞれの重みも入力として、前記パラメータを推定し、
前記（ｂ）のステップで、前記（ｄ）のステップで計算された前記残差のうち前記インライアの残差を特定し、特定した前記インライアの残差の前記統計情報に基づいて、前記閾値を計算し、
前記（ｃ）のステップで、収束していると判定する場合は、前記（ｃ）のステップにおいて、推定された前記パラメータを、前記インライアに適合するパラメータとして出力し、
前記（ｃ）のステップで、収束していないと判定する場合は、前記（ａ）のステップで推定された前記パラメータを前記初期値として、前記（ｄ）のステップ、前記（ｅ）のステップ、前記（ａ）のステップ、及び前記（ｂ）のステップそれぞれを再度実行する、
付記１１に記載のプログラム。

【0097】

（付記１３）
前記（ａ）のステップにおいて、反復再重み付け最小二乗法、又はM推定を実行して、前記パラメータを推定する、
付記１２に記載のプログラム。

【0098】

（付記１４）
前記（ｂ）のステップにおいて、前記インライアの残差の事前確率分布から前記統計情報を求め、求めた前記統計情報に基づいて、前記閾値を計算する、
付記１２または１３に記載のプログラム。

【0099】

（付記１５）
前記事前確率分布は、混合正規分布で表され、
前記（ｂ）のステップにおいて、前記統計情報として、前記インライアの残差の平均値と標準偏差とを求める、
付記１４に記載のプログラム。

【0100】

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

【産業上の利用可能性】

【0101】

以上のように本発明によれば、複数のデータ点に適合する幾何パラメータの推定において、アウトライアとインライアとを分離する際の閾値を適切に決定することができる。本発明は、複数のデータ点に適合する幾何パラメータを推定することが求められる分野、例えば、統計学、信号処理、画像処理、コンピュータビジョン等に有用である。

【符号の説明】

【0102】

１０ カメラパラメータ推定装置
１１ 残差計算部
１２ 重み更新部
１３ パラメータ推定部
１４ 閾値決定部
１５ 収束判定部
１１０ コンピュータ
１１１ ＣＰＵ
１１２ メインメモリ
１１３ 記憶装置
１１４ 入力インターフェイス
１１５ 表示コントローラ
１１６ データリーダ／ライタ
１１７ 通信インターフェイス
１１８ 入力機器
１１９ ディスプレイ装置
１２０ 記録媒体
１２１ バス

【図1】

【図2】

【図3】

【図4】

【図5】