(57)【要約】【課題】射影変換行列の最小二乗解を高速かつ精度よく推定可能とする。【解決手段】Ｍ次元の射影変換の変換行列を、該射影変換による変換前後の座標点ペアに基づいて立式される連立方程式を解くことにより導出する演算装置において、Ｎ個の座標点ペアに基づいて立式される２Ｎ個の連立方程式であって、係数行列Ａ、解ベクトルｂ、前記変換行列の各要素により構成される変数ベクトルｘによりＡｘ＝ｂと表現される連立方程式を決定する決定手段と、前記係数行列Ａを、該係数行列Ａにおける列ベクトル間の直交化を部分的に行った部分直交行列Ｑｐと上三角行列Ｒｐとに分解する分解手段と、前記部分直交行列Ｑｐの転置行列ＱｐＴを前記連立方程式の両辺に乗算して得られる変形連立方程式を処理することにより前記変数ベクトルｘを決定し、前記射影変換の変換行列の各要素を導出する導出手段と、を有する。【選択図】図４
2017059048.tif 000002
【特許請求の範囲】
【請求項１】
Ｍ次元の射影変換の変換行列を、該射影変換による変換前後の座標点ペアに基づいて立式される連立方程式を解くことにより導出する演算装置であって、
Ｎ個の座標点ペアに基づいて立式される２Ｎ個（ただし、２Ｎ＞（Ｍ＋１）２−１）の連立方程式であって、係数行列Ａ、解ベクトルｂ、前記変換行列の各要素により構成される変数ベクトルｘによりＡｘ＝ｂと表現される連立方程式を決定する決定手段と、
前記係数行列Ａを、該係数行列Ａにおける列ベクトル間の直交化を部分的に行った部分直交行列Ｑｐと上三角行列Ｒｐとに分解する分解手段と、
前記部分直交行列Ｑｐの転置行列ＱｐＴを前記連立方程式の両辺に乗算して得られる変形連立方程式を処理することにより前記変数ベクトルｘを決定し、前記射影変換の変換行列の各要素を導出する導出手段と、
を有し、
前記部分直交行列Ｑｐにおける少なくとも１組の列ベクトルは互いに直交していない
ことを特徴とする演算装置。
【請求項２】
前記分解手段は、前記係数行列Ａにおける所定の列ベクトルにおいて、それぞれの列ベクトルの有意要素の平均値を該有意要素から減算して行うことにより直交化を行う
ことを特徴とする請求項１に記載の演算装置。
【請求項３】
前記所定の列ベクトルは、前記Ｎ個の座標点ペアそれぞれに含まれる変換前の座標点の値を含む２Ｍ個の列ベクトルを含む
ことを特徴とする請求項２に記載の演算装置。
【請求項４】
前記導出手段は、前記変形連立方程式の左辺に存在するＱｐＴＱｐの非対角成分がゼロとなるように、前記連立方程式の両辺に所定の演算を行う
ことを特徴とする請求項１乃至３の何れか１項に記載の演算装置。
【請求項５】
前記導出手段は、ＱｐＴＱｐの逆行列を前記変形連立方程式の両辺に乗算し、更に、Ｒｐの逆行列を両辺に乗算することにより前記変数ベクトルｘを決定する
ことを特徴とする請求項４に記載の演算装置。
【請求項６】
前記乗算は、後退代入処理を行うことにより実現される
ことを特徴とする請求項５に記載の演算装置。
【請求項７】
Ｍ次元の射影変換の変換行列を、該射影変換による変換前後の座標点ペアに基づいて立式される連立方程式を解くことにより導出する演算装置であって、
Ｎ個の座標点ペアに基づいて立式される２Ｎ個（ただし、２Ｎ＞（Ｍ＋１）２−１）の連立方程式であって、Ｎ行（（Ｍ＋１）２−１）列の係数行列Ａｈ、Ｎ行Ｍ列の解行列Ｂ、前記変換行列の各要素により構成される（（Ｍ＋１）２−１）行Ｍ列の変数行列ＸによりＡｈＸ＝Ｂと表現される連立方程式を決定する決定手段と、
前記係数行列Ａｈを、該係数行列Ａｈにおける列ベクトル間の直交化を部分的に行った部分直交行列Ｑｐと上三角行列Ｒｐとに分解する分解手段と、
前記部分直交行列Ｑｐの転置行列ＱｐＴを前記連立方程式の両辺に乗算して得られる変形連立方程式を処理することにより前記解行列Ｂを決定し、前記射影変換の変換行列の各要素を導出する導出手段と、
を有し、
前記部分直交行列Ｑｐにおける少なくとも１組の列ベクトルは互いに直交していない
ことを特徴とする演算装置。
【請求項８】
Ｍ次元の射影変換の変換行列を、該射影変換による変換前後の座標点ペアに基づいて立式される連立方程式を解くことにより導出する演算装置における演算方法であって、
Ｎ個の座標点ペアに基づいて立式される２Ｎ個（ただし、２Ｎ＞（Ｍ＋１）２−１）の連立方程式であって、係数行列Ａ、解ベクトルｂ、前記変換行列の各要素により構成される変数ベクトルｘによりＡｘ＝ｂと表現される連立方程式を決定する決定工程と、
前記係数行列Ａを、該係数行列Ａにおける列ベクトル間の直交化を部分的に行った部分直交行列Ｑｐと上三角行列Ｒｐとに分解する分解工程と、
前記部分直交行列Ｑｐの転置行列ＱｐＴを前記連立方程式の両辺に乗算して得られる変形連立方程式を処理することにより前記変数ベクトルｘを決定し、前記射影変換の変換行列の各要素を導出する導出工程と、
を含み、
前記部分直交行列Ｑｐにおける少なくとも１組の列ベクトルは互いに直交していない
ことを特徴とする演算方法。
【請求項９】
Ｍ次元の射影変換の変換行列を、該射影変換による変換前後の座標点ペアに基づいて立式される連立方程式を解くことにより導出する演算装置における演算方法であって、
Ｎ個の座標点ペアに基づいて立式される２Ｎ個（ただし、２Ｎ＞（Ｍ＋１）２−１）の連立方程式であって、Ｎ行（（Ｍ＋１）２−１）列の係数行列Ａｈ、Ｎ行Ｍ列の解行列Ｂ、前記変換行列の各要素により構成される（（Ｍ＋１）２−１）行Ｍ列の変数行列ＸによりＡｈＸ＝Ｂと表現される連立方程式を決定する決定工程と、
前記係数行列Ａｈを、該係数行列Ａｈにおける列ベクトル間の直交化を部分的に行った部分直交行列Ｑｐと上三角行列Ｒｐとに分解する分解工程と、
前記部分直交行列Ｑｐの転置行列ＱｐＴを前記連立方程式の両辺に乗算して得られる変形連立方程式を処理することにより前記解行列Ｂを決定し、前記射影変換の変換行列の各要素を導出する導出工程と、
を有し、
前記部分直交行列Ｑｐにおける少なくとも１組の列ベクトルは互いに直交していない
ことを特徴とする演算装置。
【請求項１０】
コンピュータを、請求項１乃至７の何れか１項に記載の演算装置の各手段として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、２枚の画像間の射影変換行列を推定する技術に関するものである。
【背景技術】
【０００２】
自由度”６”を有するアフィン変換行列は、２次元画像を、拡大・縮小・回転・並行移動することが可能である。これにさらに２つの自由度を追加して水平方向と垂直方向の奥行き（遠近感）を表現できるようしたのが、２次元射影変換行列になる。変換対象を３次元空間に拡張したアフィン変換行列は自由度”１２”、射影変換行列は自由度”１５”であり、これらを使った変換処理がコンピュータグラフィックス（ＣＧ）生成で利用されている。
【０００３】
一方、コンピュータビジョンの分野では、複数の画像から変換行列を推定することも行われる。例えば、動画撮影カメラの手振れ補正を目的として、カメラ等で撮影した複数の画像から変換行列を推定し補正を行う。２つの画像中の対応点の座標データから２次元射影変換行列を求める機能は多くのアプリケーションで使われており、この変換行列はホモグラフィ行列と呼ばれる。ホモグラフィ行列は、パノラマ画像の生成、プロジェクタの投影位置の補正、平面パターンを用いるカメラ校正、対象物の形状推定・復元、道路面上の障害物検出など幅広い用途で利用することが出来る。
【０００４】
ただし、射影変換に基づく連立方程式を解いてホモグラフィ行列を求める処理は処理量が多い。そこで、特許文献１では、画像間の相関が高い場合に有効な近似演算で求める局所ホモグラフィ行列と、連立方程式を解いて求める大域的ホモグラフィ行列の２種類を使い分ける技術が開示されている。具体的には、画像間の相関が高い場合には高速演算が可能な局所ホモグラフィ行列を計算し、相関が低い場合に限り連立方程式を解いて大域的ホモグラフィ行列を求めるようにしている。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開２０１２−８６２８５号公報
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、連立方程式を作る際の対応点が多いと、係数行列Ａの行列積（ＡＴＡ）を演算する際の積和演算において、数値計算誤差が蓄積しやすくなる。一方誤差の蓄積を避ける手法として、行列のＱＲ分解を用いる方法があるが、係数行列Ａの列ベクトルを相互に直交化する必要があり非常に重い処理となる。
【０００７】
本発明は、このような問題に鑑みてなされたものであり、射影変換行列の最小二乗解を高速かつ精度よく推定可能とする技術を提供することを目的とする。
【課題を解決するための手段】
【０００８】
上述の問題点を解決するため、本発明に係る演算装置は以下の構成を備える。すなわち、Ｍ次元の射影変換の変換行列を、該射影変換による変換前後の座標点ペアに基づいて立式される連立方程式を解くことにより導出する演算装置において、Ｎ個の座標点ペアに基づいて立式される２Ｎ個（ただし、２Ｎ＞（Ｍ＋１）２−１）の連立方程式であって、係数行列Ａ、解ベクトルｂ、前記変換行列の各要素により構成される変数ベクトルｘによりＡｘ＝ｂと表現される連立方程式を決定する決定手段と、前記係数行列Ａを、該係数行列Ａにおける列ベクトル間の直交化を部分的に行った部分直交行列Ｑｐと上三角行列Ｒｐとに分解する分解手段と、前記部分直交行列Ｑｐの転置行列ＱｐＴを前記連立方程式の両辺に乗算して得られる変形連立方程式を処理することにより前記変数ベクトルｘを決定し、前記射影変換の変換行列の各要素を導出する導出手段と、を有し、前記部分直交行列Ｑｐにおける少なくとも１組の列ベクトルは互いに直交していない。
【０００９】
又は、Ｍ次元の射影変換の変換行列を、該射影変換による変換前後の座標点ペアに基づいて立式される連立方程式を解くことにより導出する演算装置において、Ｎ個の座標点ペアに基づいて立式される２Ｎ個（ただし、２Ｎ＞（Ｍ＋１）２−１）の連立方程式であって、Ｎ行（（Ｍ＋１）２−１）列の係数行列Ａｈ、Ｎ行Ｍ列の解行列Ｂ、前記変換行列の各要素により構成される（（Ｍ＋１）２−１）行Ｍ列の変数行列ＸによりＡｈＸ＝Ｂと表現される連立方程式を決定する決定手段と、前記係数行列Ａｈを、該係数行列Ａｈにおける列ベクトル間の直交化を部分的に行った部分直交行列Ｑｐと上三角行列Ｒｐとに分解する分解手段と、前記部分直交行列Ｑｐの転置行列ＱｐＴを前記連立方程式の両辺に乗算して得られる変形連立方程式を処理することにより前記解行列Ｂを決定し、前記射影変換の変換行列の各要素を導出する導出手段と、を有し、前記部分直交行列Ｑｐにおける少なくとも１組の列ベクトルは互いに直交していない。
【発明の効果】
【００１０】
本発明によれば、射影変換行列の最小二乗解を高速かつ精度よく推定可能とする技術を提供することができる。
【図面の簡単な説明】
【００１１】
【図１】ホモグラフィ行列を求めるための連立方程式の係数行列と転置行列を例示的に示す図である。
【図２】正規方程式の係数行列を例示的に示す図である。
【図３】第１実施形態における係数行列の変遷を説明する図である。
【図４】第１実施形態における処理のフローチャートである。
【図５】第２実施形態における係数行列の変遷を説明する図である。
【図６】第２実施形態における処理のフローチャートである。
【図７】第３実施形態における係数行列の変遷を説明する図である。
【図８】第３実施形態における処理のフローチャートである。
【図９】ＰＣにより構成した演算装置の内部ブロック図である。
【発明を実施するための形態】
【００１２】
以下に、図面を参照して、この発明の好適な実施の形態を詳しく説明する。なお、以下の実施の形態はあくまで例示であり、本発明の範囲を限定する趣旨のものではない。
【００１３】
（第１実施形態）
本発明に係る演算装置の第１実施形態として、ホモグラフィ行列や３次元の射影変換行列等の最小二乗解を高速かつ精度よく推定する演算装置について説明する。具体的には、連立方程式の係数行列の構造に着目し、列ベクトルの直交化を部分的に行って最小二乗推定を行う。これにより高精度かつ高速な推定を可能とする。
【００１４】
＜前提技術＞
まず、ホモグラフィ行列を利用して実現する動画撮影カメラの手振れ補正の原理、及び、ホモグラフィ行列の求め方を説明する。
【００１５】
動画撮影して得られる連続フレームのフレーム間のブレを幾何変換行列として検出し、この幾何変換行列を十数フレームに渡って平滑化した行列を導出する。このようにして導出された行列を用いて画像を逆補正することにより、フレーム間のブレや視点移動を平滑化でき、動画撮影時の手ぶれ補正機能を実現できる。特に、カメラから被写体までの距離が十分に離れており、被写体の各注目点が同一平面上に存在すると見なすことが出来る場合、フレーム間のブレを表現する幾何変換行列はホモグラフィ行列で表すことができる。そのとき、時間的に先行するフレームと後続のフレームとの間の画素の座標は式（１）で対応付けられる。
【００１６】
【数１】
2017059048.tif 000003
【００１７】
式（１）における右辺の３行３列の行列がホモグラフィ行列であり、先行するフレームの（ｕ，ｖ）座標の画素（変換前の座標点）は、（ｘ，ｙ）座標の画素（変換後の座標点）に変換される。以下では、これらの変換前後の２つの座標点を座標点ペアと呼ぶ。また、式（１）は、非線形な射影変換を線形演算らしく表現するために、同次座標を用いた変換式で表現される。
【００１８】
右下の９番目の要素が”１”となるように規格化（正規化）すると、他の８つの要素が自由な値をとることが可能であり、自由度が８の行列となる。一方、式（１）の関係で表される２画像間の座標を展開すると、式（２）に示す２つの方程式が得られる。
【００１９】
【数２】
2017059048.tif 000004
【００２０】
４組の座標点ペアから８つの方程式、すなわち、８元連立方程式を決定し立式することができ、この連立方程式を解くことによりホモグラフィ行列を１つ求めることができる。但し、４組の座標点ペアは被写体空間において同一平面上に存在する必要があり、同一平面上に無い座標点ペアを含む場合、意味の無いホモグラフィ行列になってしまう。そのため、どれだけ適切な解であるかを、他の座標点ペアの座標データを変換して評価する必要がある。
【００２１】
一方、５組以上の座標点からホモグラフィ行列を求める場合、方程式の数が８つを超える過条件の連立方程式となり、全ての方程式を満たす解が存在しなくなり、方程式を解く事ができなくなってしまうことがある。Ｍ次元の射影変換においては、Ｎ個以上（ただし、２Ｎ＞（Ｍ＋１）２−１）の座標点ペアを用いる場合に２Ｎ個の連立方程式が立式され過条件の連立方程式となる。
【００２２】
この場合は、方程式の係数行列Ａを転置した行列ＡＴを、両辺に、すなわち係数行列Ａと解ベクトルｂに乗じる。これにより、係数行列と解ベクトルを８次元空間に投影し、「正規方程式」と呼ばれる８元連立方程式に変換して解を算出する。変数ベクトルをｘで表した時、正規方程式は（ＡＴＡ）ｘ＝ＡＴｂと表現される。ここで（ＡＴＡ）はの正規方程式の係数行列であり、行列のサイズは８×８である。右辺のＡＴｂは８次元ベクトルである。
【００２３】
正規方程式を解いて得られる解は、５組以上の座標点ペアに対する最小二乗法解（最小二乗法に基づく解）であり、全座標点ペアを反映した解となる。この正規方程式を用いて最小二乗解を得る手法は、統計処理でよく用いられる一般的な手法である。
【００２４】
＜装置構成＞
図９は、演算装置のハードウェア構成を示すブロック図である。ここでは、演算装置１０１が一般的なパーソナルコンピュータ（ＰＣ）により構成する例を示しているためディスプレイ１２３、通信Ｉ／Ｆ１２４、キーボード１２５などを含んでいるが、これらは必須の構成というわけではない。
【００２５】
ＣＰＵ１２０は、各種演算を実行するほか、演算装置１０１を統括的に制御する。ＣＰＵ１２０は、例えばＲＯＭ１２２やハードディスクドライブ（ＨＤＤ）１２６に格納された演算プログラムを実行することにより図４に示される各処理を実現する。
【００２６】
ＨＤＤ１２６、例えば、演算装置１０１で利用される演算プログラムあるいは各種の制御プログラムを記憶する。また、演算プログラムあるいは各種の制御プログラムに関連する各種情報を保存する。ここでは、各種情報としては、連立方程式に関連する情報を含む。また、一時的に各種情報を記憶するためにＲＡＭ１２１も利用される。
【００２７】
キーボード１２５は、ユーザから、計算に用いる各種情報などの入力を受け付ける機能部である。また、ディスプレイ１２３は、ユーザに対して、計算処理の進行状況や計算結果など各種情報を提供する機能部である。通信インタフェース（Ｉ／Ｆ）１２４は、不図示の外部装置に接続するためのインタフェースであり、例えば、有線通信又は無線通信のインタフェースである。Ｉ／Ｆ１２４を介して、例えば、外部の撮像装置などから処理の対象となる複数フレームの画像を入力したり、特徴点のペアを入力したりすることができる。また、Ｉ／Ｆ１２４を介して、得られた射影変換行列を外部の画像補正装置などに出力することができる。
【００２８】
＜ホモグラフィ行列導出の動作概要＞
まず、Ｎ組（Ｎは５以上）の座標点ペアに対する過条件連立方程式と、その正規方程式の具体例を示す。既に説明したように、Ｎ組の座標点ペアから２Ｎ個の８元連立方程式を決定することが出来る。この連立方程式は、係数行列Ａ、変数ベクトルｘ、解ベクトルｂを用いてＡｘ＝ｂと表現される。
【００２９】
図１は、ホモグラフィ行列を求めるための連立方程式の係数行列・解ベクトルと、係数行列の転置行列を例示的に示す図である。図１（Ａ）は、２Ｎ個の８元連立方程式を、行列形式で表現したものである。図１（Ｂ）は、図１（Ａ）の行列を解ベクトルに基づき並び替えたものである。以下では、図１（Ｂ）に示す行列および列ベクトルをそれぞれ係数行列Ａ、解ベクトルｂと呼ぶ。また、図１（Ｃ）は、係数行列Ａの転置行列であるＡＴを示している。
【００３０】
図２は、正規方程式の係数行列を例示的に示す図である。図２（Ａ）に示す行列と図２（Ｂ）に示す列ベクトルは、係数行列Ａと解ベクトルｂの各々に、図１（Ｃ）の転置行列を掛けて８次元空間へ投影したものである。
【００３１】
すなわち、各々は正規方程式の両辺に存在する係数行列（ＡＴＡ）と解ベクトル（ＡＴｂ）に対応する。Σ記号は不図示の添え字ｉ＝１，２，・・・，Ｎの合算を表す。なお、この係数行列は対称行列になる。つまり、図２（Ａ）においては、左下の２行６列をΣ記号のみを用いて記述しているが、右上の２列６行を転置した要素と同じである。
【００３２】
図２（Ａ）を見ると、ホモグラフィ行列を求める正規方程式にはゼロ要素が規則的に配置されており、その係数行列には冗長な部分が存在していることが分かる。
【００３３】
ＱＲ分解法では係数行列Ａを、直交行列Ｑと上三角行列Ｒとに分解する。一方、第１実施形態では、係数行列Ａを、部分的に直交させた部分直交行列ＱｐとＱｐに対応する上三角行列Ｒｐとに分解する。そして、Ｑｐの転置行列をＱｐに乗算して、得られる８×８行列のゼロ要素を、図２（Ａ）の行列に比較してさらに増やす。これにより、図２に示す正規方程式より解き易い形の方程式にする。
【００３４】
特に、第１実施形態では、図１（Ｂ）に示す係数行列中の１列目と４列目を除く４個の列ベクトルにおいて、上側Ｎ個と下側Ｎ個に対し、各々の平均を算出して、各要素からこれら平均を減算する。“０”がＮ個並んでいるところは何もしなくてもよい。なお、Ｍ次元の射影変換においては、Ｎ個の座標点ペアそれぞれに含まれる変換前の座標点の値を含む２Ｍ個の列ベクトルに対して同様の処理を行うことになる。
【００３５】
平均値を減算した各列ベクトルは、１列目及び４列目と直交するようになる。このように処理した行列をＱｐとする。すなわち、一般的なＱＲ分解におけるＱ行列は、列ベクトルが相互に直交している。それに対し、第１実施形態におけるＱｐの列ベクトルは特定の列に対して直交しているだけであり、互いに直交していない列ベクトルの組が存在する。
【００３６】
図３は、第１実施形態における係数行列の変遷を説明する図である。具体的には、連立方程式である正規方程式の両辺に対して同様の処理を行い、変形連立方程式を生成する処理に相当する。図３（Ａ）は、係数行列Ａを、上述のＱｐと上三角行列Ｒｐの積で表したものである。すなわち、図３（Ａ）の右側に示す上三角行列がＲｐである。なお、図３においては、ｕ座標、ｖ座標から上述の平均を減算した結果をｄｕ、ｄｖと表記している。
【００３７】
次に、Ｑｐの転置行列ＱｐＴを左辺のＱｐと右辺の列ベクトルｂに乗算する。このとき、ＱｐＴＱｐの計算結果は８×８の対称行列で、図３（Ｂ）に示すような行列となる。対称行列であるため、左下三角行列部分の記述を省略した。そして、この８×８の係数行列に対して、以下の処理を行う。なお以下の処理は、結果的には、ＱｐＴＱｐの逆行列を乗算する処理に相当する。
【００３８】
まず、Ｎの逆数、すなわち１／Ｎを算出して、１行目と４行目に乗算する。次に、２，３行目と５，６行目の対角要素近辺の２×２領域を対角化するため、２×２行列の逆行列を求め、この逆行列を用いて、２，３行目と５，６行目の各列を２次の列ベクトルとして変換する。なお、対角要素近辺の２×２領域は、変換すると単位行列になることが分かっているため、変換は不要である。
【００３９】
すなわち、上述の逆行列により実際に変換しなければならない列ベクトルは、７列目、８列目と解ベクトルにのみ存在する。この２×２逆行列による変換によって、８×８係数行列は、左上６×６領域を単位行列化した行列となる。
【００４０】
次に、７，８行目の１〜６列の要素をゼロに消去する。７行目の４つの有意要素と、７列目の４つの有意要素との間で積和を演算し、行と列が交差する７行７列目の要素から減算する。この積和・減算処理を、７行目と８列目、８行目と７列目、８行目と８列目の間で同様に行う。
【００４１】
以上の処理により、７，８行における１〜６列の要素がゼロ消去される。次は、残りの７，８列目の２×２要素を対角化する。ここでは、上述の２×２の逆行列で、対応する解ベクトルの要素を変換するだけでよい。これにより変数ベクトルｘの各要素を決定することが可能となる。
【００４２】
＜演算装置の動作＞
図４は、第１実施形態における処理のフローチャートである。以下の処理は、例えば、時間的に先行するフレームと後続のフレームとの間の座標点ペア群が入力されることにより開始される。
【００４３】
ステップＳ４０１では、ＣＰＵ１２０は、座標点ペア群の中からＮ組（Ｎは５以上）の座標データを選択し、図１（Ｂ）に示す形式で係数行列Ａと解ベクトルｂを生成する。
【００４４】
ステップＳ４０２では、ＣＰＵ１２０は、２Ｎ行×８列の係数行列Ａを、図３（Ａ）に示す形式で行列Ｑｐと、行列Ｒｐとに分解する。上述のように、Ｑｐは、係数行列Ａの１列目と４列目を除く列ベクトルにおいて、上側Ｎ個と下側Ｎ個に対し各々の平均を算出して減算することにより得られる。
【００４５】
ステップＳ４０３では、ＣＰＵ１２０は、Ｑｐの転置行列ＱｐＴを、Ｑｐとｂに乗算する。
【００４６】
ステップＳ４０４では、ＣＰＵ１２０は、Ｎの逆数と、２×２行列の逆行列を計算する。ステップＳ４０５では、ＣＰＵ１２０は、図３（Ｂ）に示す行列に対し、１，４行目にＮの逆数を、２，３行目と５，６行目に２×２逆行列を乗算する。ステップＳ４０６では、ＣＰＵ１２０は、図３（Ｂ）に示す行列に対し、７，８行目の１〜６列をゼロ消去する。
【００４７】
ステップＳ４０７では、ＣＰＵ１２０は、７，８行目に存在する残りの係数から２元連立方程式を解き、パラメータを２つ算出する。ステップＳ４０８では、ＣＰＵ１２０は、Ｑｐの７，８列に残る要素を消去する後退代入処理を行う。ステップＳ４０９では、ＣＰＵ１２０は、上三角行列Ｒｐの非対角要素を消去する後退代入処理を行う。
【００４８】
これらの一連の処理により、最小二乗推定したホモグラフィ行列を１個得ることが出来る。なお、上述の処理では、ゼロでない有意な係数のみを効率よく処理するため、細々とした処理ステップに分割した。ただし、図３（Ｂ）に示すＱｐＴＱｐをガウス消去法やガウス・ジョルダン法などのよく知られた方法で処理することも可能である。その場合は、上述のＳ４０４〜Ｓ４０８をこれらの方法に置き換えて処理することになる。また、上述の説明では、ゼロ要素も区別なく処理しているが、ゼロ要素を検出して処理をスキップするよう構成してもよい。
【００４９】
なお、Ｓ４０１で選択した座標点ペアの中には、前提から外れたアウトライアが含まれている可能性がある。なお、アウトライアとは、先行するフレームにおける座標をホモグラフィ行列により変換した結果の座標と、対応する後続のフレームにおける座標とが所定閾値よりも大きい差分（ズレ）を有する座標点ペアを意味する。
【００５０】
そのため、一部または全部の座標データを入れ替えて、いくつもの最小二乗推定行列を得るように構成してもよい。すなわち、種々の組合せのＮ組の座標点ペアについて図４に示す処理を行い、変換によるズレ具合の評価を行い、評価の一番良いパラメータを選んで、最終的なホモグラフィ行列として導出してもよい。
【００５１】
以上説明したとおり第１実施形態によれば、直交化処理が容易な一部の列ベクトル間のみを直交化することにより、ホモグラフィ行列や３次元の射影変換行列等の最小二乗解を高速かつ精度よく推定可能となる。すなわち、少なくとも１組の列ベクトル間は直交処理を行わない。これにより、単位時間内に算出した多くのパラメータを評価して、信頼性の高い変換行列を選別することが可能となる。その結果、例えば、動画像データを処理する場合には、１フレーム時間内に多くのパラメータを評価することができるため、より精度の高いホモグラフィ行列を導出することが可能となる。その他、画像合成や位置合せ、カメラ防振、対象物の３次元形状推定・復元などの性能を向上することが可能となる。
【００５２】
（第２実施形態）
第２実施形態では、第１実施形態と同等に、係数行列Ａを、部分的に直交させた部分直交行列Ｑｐと該Ｑｐに対応する上三角行列Ｒｐとに分解する。ただし、Ｑｐの算出方法が第１実施形態と異なる。
【００５３】
具体的には、係数行列Ａにおいて、２列目，３列目に対し部分平均を算出・減算して１列目に直交化させた後、さらに２列目と３列目を直交化させて、部分直交行列Ｑｐとする。この場合、４，５，６列目はＱｐＴＱｐを計算する時、１，２，３列目の結果を利用すれば足りるため処理を省略できる。
【００５４】
図５は、第２実施形態における係数行列の変遷を説明する図である。図５（Ａ）は、係数行列Ａを、上述のＱｐと上三角行列Ｒｐの積で表したものである。なお、図３と同様に、ｕ座標、ｖ座標から上述の部分平均を減算した結果をｄｕ、ｄｖと表記している。
【００５５】
図５（Ｂ）は、ＱｐＴＱｐにより得られる行列を示している。７，８列目の非零の要素が増えているが、左上６×６の部分行列が、ＱｐＴＱｐ計算直後の段階で既に対角行列になっている。よって、図５（Ｂ）に示す行列においては、１〜６行目の各行を対角要素の値で除算することによって、この６×６部分行列を単位行列化することができる。
【００５６】
次に、図５（Ｂ）に示す行列の７，８行目の１〜６列の要素をゼロに消去する。まず、７行目１〜６列の６つの要素と、７列目１〜６行の６つの要素間で積和を演算し、行と列が交差する７行７列目の要素から減算する。この積和・減算処理を、７行目と８列目、８行目と７列目、８行目と８列目の間で同様に行う。以上の処理により、７，８行における１〜６列の要素がゼロ消去される。以降は、第１実施形態と同様の処理であるため、説明を省略する。
【００５７】
＜演算装置の動作＞
図６は、第２実施形態における処理のフローチャートである。第１実施形態と同様に、例えば、時間的に先行するフレームと後続のフレームとの間の座標点ペア群が入力されることにより開始される。なお、ステップＳ６０２、Ｓ６０４、Ｓ６０５の３つのステップ以外は第１実施形態と同様であるため説明は省略する。
【００５８】
ステップＳ６０２では、ＣＰＵ１２０は、２Ｎ行×８列の係数行列Ａを、図５（Ａ）に示す形式で行列Ｑｐと、行列Ｒｐとに分解する。上述のように、Ｑｐは、２列目，３列目に対し部分平均を算出・減算して１列目に直交化させた後、さらに２列目と３列目を直交化させることにより得られる。上述したように４，５，６列目は、１，２，３列目の結果を利用すれば足りる。
【００５９】
ステップＳ６０４では、ＣＰＵ１２０は、図５（Ｂ）に示す行列（ＱｐＴＱｐの結果）における、対角行列になっている左上６×６行列の対角要素の逆数を計算する。なお、対角要素は計６つあるが、４〜６行の対角要素は１〜３行と同じ値になるため、実際には３つの逆数のみ計算すればよい。
【００６０】
ステップＳ６０５では、ＣＰＵ１２０は、Ｓ６０４により計算された３つの逆数を１〜６行の対応する行に乗算する。実際には、行の要素すべてに乗算する必要は無く、７，８列目と解ベクトルの３つの要素に乗算するだけでよい。
【００６１】
なお、第２実施形態においては、７列目と８列目の列ベクトルの直交化を行わなかったが、第１実施形態と同様に、この２列を１列目と４列目に対して直交化しておくことも可能である。その場合の直交化以降の処理は、第１実施形態とほぼ同じであり、積和を計算するときの要素数（＝積の回数）が６つから４つに減少する点のみが異なる。
【００６２】
以上説明したとおり第２実施形態によれば、直交化処理が容易な一部の列ベクトル間のみを直交化することにより、ホモグラフィ行列や３次元の射影変換行列等の最小二乗解を高速かつ精度よく推定可能となる。
【００６３】
（第３実施形態）
第３実施形態では、より小さな行列サイズでホモグラフィ行列や３次元の射影変換行列等の最小二乗解を高速かつ精度よく推定する形態について説明する。特に、行列サイズを半分に減らした場合の係数データの保持方法と処理方法について説明する。
【００６４】
図１（Ｂ）に示すように、列ベクトルを直交化する対象である係数行列Ａは２Ｎ行×８列のサイズである。つまり、Ｎの数が例えば５００である場合には係数行列を格納するための記憶領域は８Ｋワードになり、かなり大きくなる。そのため、メモリサイズをあまり大きく取れない場合には適切に処理できない可能性がある。特に、組込み用途のプロセッサなどでは、メモリ領域を他の処理と分割して使用することが想定されるため、行列サイズは小さい程よい。
【００６５】
図７は、第３実施形態における係数行列の変遷を説明する図である。図７（Ａ）は、係数行列Ａｈと解行列Ｂを示す図である。係数行列Ａｈは、行数を図１（Ｂ）の半分のＮ行×８列とした行列であり、行列へ座標データを格納して生成する。Ｍ次元の場合は、係数行列Ａｈは、Ｎ行（（Ｍ＋１）２−１）列となる。また、図１（Ｂ）に示す解ベクトルｂのＮ＋１行目以降を２列目に移動して、Ｎ行×２列の解行列Ｂとする。Ｍ次元の場合は、解行列Ｂは、Ｎ行Ｍ列となる。また、変数行列Ｘは、（（Ｍ＋１）２−１）行Ｍ列となる。このとき、連立方程式は、ＡｈＸ＝Ｂと表現されることになる。
【００６６】
ここで、係数行列Ａｈの１列目と４列目を除く各列を１列目と直交化する。すなわち、各列のＮ個のデータの平均を算出し各データから減算・除去して部分直交行列Ｑｐとする。そして、対応する８行×８列の上三角行列Ｒｐを生成する。
【００６７】
図７（Ｂ）は、係数行列Ａｈを、上述のＱｐと上三角行列Ｒｐの積で表したものである。なお、図３と同様に、ｕ座標、ｖ座標から上述の平均を減算した結果をｄｕ、ｄｖと表記している。
【００６８】
図７（Ｃ）は、ＱｐＴＱｐにより得られる８行×８列の行列を示している。ここでは、図７（Ｃ）においてΣ記号を付した１７個の要素を計算することになる。左上隅の要素の値“Ｎ”は確定値であるため計算対象には含まれない。
【００６９】
この１７個の要素の中で５，６列目のｙを含む７つの要素を、７，８列目に移動する。もしくは、移動した先の要素に加算して図３（Ｂ）に示す行列の７，８列目と同じ構成の列ベクトルを形成する。但し、７列目の一番下の要素は無視する。その後、行列を対称化するため、該行列の上三角領域を対角成分で折り返して下三角領域へコピーする。
【００７０】
図７（Ｄ）は、ＱｐＴを解行列Ｂに乗算して得られる８行×２列の行列である。図７（Ｄ）においてΣ記号を付した１０個の要素を計算し、上述と同様に移動と加算を行って、図２（Ｂ）に示すベクトルと同じ結果を得る。
【００７１】
これ以降は、第１実施形態における処理フローのＳ４０４以降と同じ処理になる。これは、同じ図３（Ｂ）の行列を処理するためである。
【００７２】
＜演算装置の動作＞
図８は、第３実施形態における処理のフローチャートである。上述のように、Ｓ４０４以降は第１実施形態と同じであるため、それ以前のステップＳ８０１〜Ｓ８０５について説明する。
【００７３】
ステップＳ８０１では、ＣＰＵ１２０は、座標点ペア群の中からＮ組（Ｎは５以上）の座標データを選択し、図７（Ａ）に示す形式でＮ行×８列の係数行列Ａｈを生成する。
【００７４】
ステップＳ８０２では、ＣＰＵ１２０は、係数行列Ａｈを、図７（Ｂ）に示す形式で行列Ｑｐと、行列Ｒｐとに分解する。上述のように、Ｑｐは、係数行列Ａｈの１列目と４列目を除く各列から列毎の平均値を減算することに得られる。
【００７５】
ステップＳ８０３では、ＣＰＵ１２０は、Ｑｐの転置行列ＱｐＴを、Ｑｐに乗算する。上述したように、形成される行列は８行×８列の行列であるが、この段階では１７個の要素のみを計算する。
【００７６】
ステップＳ８０４では、ＣＰＵ１２０は、形成中の８×８行列を、図３（Ｂ）に示す行列に類似するように変形する。例えば、上述したように、５，６列目のｙを含む７つの要素を、７，８列目に移動するか、もしくは、移動した先の要素に加算する。ステップＳ８０５では、ＣＰＵ１２０は、行列の上三角領域を対角成分で折り返して下三角領域へコピーし、行列を対称化する。
【００７７】
以上の処理の結果、図３（Ｂ）に示す８行×８列の対称行列が得られることになる。そのため、この処理結果に対してＳ４０４〜Ｓ４０９の処理を行えば、第１実施形態と同様の結果を得ることが出来る。
【００７８】
以上説明したとおり第３実施形態によれば、係数行列を格納する行列のサイズを小さくして記憶領域を減らすことが可能となる。なお、次のようにすれば、行列のサイズをさらに小さくすることができる。
【００７９】
係数行列Ａｈにおいては、４列目に格納されるデータは全部“０”であるため、該４列目を除いて７列にすることも考えられる。また、１列目に格納するデータは全部“１”であるため、該１列目を除いて６列にすることも考えられる。
【００８０】
その場合、転置行列との乗算によって得られる行列サイズは６行×６列となる。図３（Ｂ）に示す８行×８列の行列を得るには、要素を大幅に並び替える必要があるものの、係数行列を格納するための記憶領域をさらに２５％減らすことが可能となる。
【００８１】
（変形例）
上述の第１〜第３実施形態では、２次元の射影変換行列であるホモグラフィ行列を対象として説明を行った。ただし、本発明は、Ｍ次元、例えば３次元の射影変換行列の最小二乗推定にも適用することが可能である。
【００８２】
ホモグラフィ行列の自由度は８（＝３×３−１）であるため、Ｎ個の座標点ペアに対する係数行列Ａのサイズは２Ｎ行×８列である。そして、ＱｐＴＱｐの演算結果は８行×８列の行列である。
【００８３】
一方、３次元射影変換行列の自由度は１５（＝４×４−１）であるため、係数行列のサイズは２Ｎ行×１５列となる。つまり、対応するＱｐＴＱｐの演算結果は１５行×１５列の行列になる。そのため、１５列の係数行列の全ての列ベクトルをＱＲ分解を用いて完全に直交化する場合、８列の場合の４倍近い計算量となる。
【００８４】
一方、上述の実施形態と同様に、平均値除去を中心とする部分直交化を行った場合は、短時間で処理することが可能である。すなわち、ＱｐＴＱｐの結果は１５×１５行列であるが、非ゼロの非対角成分を効率的にゼロ消去することが可能である。
【００８５】
（その他の実施例）
本発明は、上述の実施形態の１以上の機能を実現するプログラムを、ネットワーク又は記憶媒体を介してシステム又は装置に供給し、そのシステム又は装置のコンピュータにおける１つ以上のプロセッサーがプログラムを読出し実行する処理でも実現可能である。また、１以上の機能を実現する回路（例えば、ＡＳＩＣ）によっても実現可能である。
【符号の説明】
【００８６】
１０１  演算装置；  １２０  ＣＰＵ；  １２１  ＲＡＭ；  １２２  ＲＯＭ
【図１】
2017059048.tif 000005
【図２】
2017059048.tif 000006
【図３】
2017059048.tif 000007
【図４】
2017059048.tif 000008
【図５】
2017059048.tif 000009
【図６】
2017059048.tif 000010
【図７】
2017059048.tif 000011
【図８】
2017059048.tif 000012
【図９】
2017059048.tif 000013