(57)【要約】【課題】演算負荷量を削減して、整数グリッド選択を柔軟にし、かつ、サブピクセル精度の移動位置推定精度を高めるようにした画像処理装置を提供する。【解決手段】画像処理装置の画像入力手段は、複数の画像を入力し、変換手段は、前記画像入力手段によって入力された画像を整数グリッド上で変換し、係数算出手段は、相違度又は類似度を求める関数を近似する２次曲面の係数を算出し、推定手段は、前記係数算出手段によって算出された係数のサブ量子化精度での推定を行う。【選択図】図１
2008209971.tif 000002
【特許請求の範囲】
【請求項１】
複数の画像を入力する画像入力手段と、
前記画像入力手段によって入力された画像を整数グリッド上で変換する変換手段と、
相違度又は類似度を求める関数を近似する２次曲面の係数を算出する係数算出手段と、
前記係数算出手段によって算出された係数のサブ量子化精度での推定を行う推定手段
を具備することを特徴とする画像処理装置。
【請求項２】
前記係数算出手段は、近似の場合に各整数グリッド位置に関する重み付けを行う
ことを特徴とする請求項１に記載の画像処理装置。
【請求項３】
前記係数算出手段による重み付けは、各整数グリッドを中心とするボロノイ分割領域の面積又は超体積である
ことを特徴とする請求項２に記載の画像処理装置。
【請求項４】
前記ボロノイ分割を行う領域は、２次曲面で近似を行いたい領域内のみに限定する
ことを特徴とする請求項３に記載の画像処理装置。
【請求項５】
（５２）式と同等のものが存在するか否かを判断する判断手段
【数５２】
2008209971.tif 000091
をさらに具備し、
前記判断手段によって存在しないと判断された場合には、前記ボロノイ分割を行う領域を広げる
ことを特徴とする請求項４に記載の画像処理装置。
【請求項６】
前記領域を広げる場合は、全軸で対称となるように１整数グリッドずつ広げる
ことを特徴とする請求項５に記載の画像処理装置。
【請求項７】
整数グリッド上で、整数の範囲で相違度が最小となる点を探索する探索手段
をさらに具備することを特徴とする請求項６に記載の画像処理装置。
【請求項８】
前記変換手段による変換は、Ｎ次元であり、Ｎ変数の２次曲面の係数を推定し、
さらに、
（２５）式における行列Ａと同等のものが半正定値であるか否かを検証する検証手段と、
【数２５】
2008209971.tif 000092
Ｍｏｏｒ−Ｐｅｎｒｏｓｅの一般化逆行列と同等のものを算出する一般化逆行列算出手段と、
サブ量子化精度での変換係数を算出する変換係数算出手段
を具備することを特徴とする請求項１に記載の画像処理装置。
【請求項９】
前記変換手段による変換は、Ｎ次元であり、Ｎ変数の２次曲面の係数を推定し、
さらに、
（２５）式における行列Ａと同等のものが正定値であるか否かを検証する検証手段と、
【数２５】
2008209971.tif 000093
サブ量子化精度での変換係数を算出する変換係数算出手段
を具備することを特徴とする請求項１に記載の画像処理装置。
【請求項１０】
前記変換手段による変換は、平行移動であり、２変数の２次曲面の係数を推定し、
さらに、
（８）式と同等のものが成り立つか否かを検証する検証手段と、
【数８】
2008209971.tif 000094
（９）式と同等のものを算出する算出手段
【数９】
2008209971.tif 000095
を具備することを特徴とする請求項１に記載の画像処理装置。
【請求項１１】
さらに、
（１０）式と同等のものを検証する検証手段と、
【数１０】
2008209971.tif 000096
（１２）式と同等のものを算出する算出手段
【数１２】
2008209971.tif 000097
を具備することを特徴とする請求項８、９又は１０に記載の画像処理装置。
【請求項１２】
相違度又は類似度を算出する整数グリッド位置を固定として、予め（５３）式と同等のもの又は（５４）式と同等のものを計算し、（５５）式と同等のもの又は（５６）式と同等のもので２次曲面の係数を推定する
【数５３】
2008209971.tif 000098
【数５４】
2008209971.tif 000099
【数５５】
2008209971.tif 000100
【数５６】
2008209971.tif 000101
ことを特徴とする請求項１に記載の画像処理装置。
【請求項１３】
相違度又は類似度を算出する整数グリッド位置は、縦３×横３の９点とし、（３５）式と同等のものと等価の式を用いて、２次曲面の係数を推定する
【数３５】
2008209971.tif 000102
ことを特徴とする請求項１に記載の画像処理装置。
【請求項１４】
相違度又は類似度を算出する整数グリッド位置は、縦３×横３の９点とし、（４０）式と同等のものと等価の式を用いて、２次曲面の係数を推定する
【数４０】
2008209971.tif 000103
ことを特徴とする請求項１に記載の画像処理装置。
【請求項１５】
さらに、
相違度又は類似度を算出した整数グリッド位置を入力するグリッド位置入力手段と、
行列Ｂを算出する行列算出手段と、
連立一次方程式（１９）式と同等のものを計算する計算手段
【数１９】
2008209971.tif 000104
を具備することを特徴とする請求項１２、１３又は１４に記載の画像処理装置。
【請求項１６】
さらに、
行列Ｂと同等のもののランクを算出する又は（７１）式と同等のものの正則性を検証する検証手段と、
【数７１】
2008209971.tif 000105
相違度を算出する整数グリッドを増加させる整数グリッド増加手段
を具備することを特徴とする請求項１５に記載の画像処理装置。
【請求項１７】
さらに、
整数グリッド位置と、行列（５３）式と同等のもの又は行列（５４）式と同等のものの関係を保持する関係保持手段と、
入力された整数グリッド位置が前記関係保持手段に保持されているか否かを判断する判断手段
【数５３】
2008209971.tif 000106
【数５４】
2008209971.tif 000107
を具備し、
前記係数算出手段は、テーブル内の行列（５３）式と同等のもの又は行列（５４）式と同等のものを用いて、２次曲面の係数を算出する
ことを特徴とする請求項１に記載の画像処理装置。
【請求項１８】
前記判断手段によって、入力された整数グリッド位置が前記関係保持手段に保持されていないと判断された場合、整数グリッドを増加させる整数グリッド増加手段
をさらに具備することを特徴とする請求項１７に記載の画像処理装置。
【請求項１９】
前記推定手段は、画像ブロックの整数移動位置の相違度又は類似度を複数算出し、それを用いてサブピクセル精度の動きベクトルを推定する
ことを特徴とする請求項１に記載の画像処理装置。
【請求項２０】
さらに、
全ての変換係数を０．５グリッド分ずらす画像を作成する画像作成手段
を具備し、
前記推定手段は、前記画像作成手段によって生成された画像で求めた位置と、ずらす前の画像で求めた位置とから、（４５）式と同等のもので位置を求める
【数４５】
2008209971.tif 000108
ことを特徴とする請求項１に記載の画像処理装置。
【請求項２１】
さらに、
一つ以上の変換係数を０．５グリッド分ずらす画像を作成する画像作成手段
を具備し、
前記推定手段は、前記画像作成手段によって生成された画像で求めた位置と、ずらす前の画像で求めた位置とから、（４７）式と同等のもの、（４８）式と同等のもの又は（４９）式と同等のものを用いた式で位置を求める
【数４７】
2008209971.tif 000109
【数４８】
2008209971.tif 000110
【数４９】
2008209971.tif 000111
ことを特徴とする請求項１に記載の画像処理装置。
【請求項２２】
画像処理システムを、
複数の画像を入力する画像入力手段と、
前記画像入力手段によって入力された画像を整数グリッド上で変換する変換手段と、
相違度又は類似度を求める関数を近似する２次曲面の係数を算出する係数算出手段と、
前記係数算出手段によって算出された係数のサブ量子化精度での推定を行う推定手段
として機能させることを特徴とする画像処理プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、画像処理装置及び画像処理プログラムに関する。
【背景技術】
【０００２】
複数の画像の位置合わせを行うこと、つまり複数の画像内で対応する一部の画像を重ね合わせることによって１枚の画像を生成することにより様々な応用が考えられる。
例えば、次の３つのような応用がある。
（１）画像を貼り合わせてモザイク画像を作成する。
（２）ステレオ画像のマッチングを取ることによって、３次元空間認識を行う。
（３）動画像の動き補償ベクトルを求めることによって、動画像圧縮率を高めることができる。
特に、位置合わせの精度をサブピクセル精度とすることによって、モザイク画像の画質、３次元空間認識の精度、あるいは、動画像圧縮率をさらに高めることができる。
【０００３】
例えば、特許文献１には、画像の位置合わせをサブピクセル精度で行う技術が記載されている。すなわち、Ｎ次元類似度を考慮することで、画像間対応パラメータを少ない計算量で安定かつ高速に高精度で同時に推定できるようにした、画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法及び多パラメータ高精度同時推定プログラムを提供することを課題とし、離散的な位置で得られた画像間のＮ次元類似度値が、ある１つのパラメータ軸に対して平行な直線上で最大又は最小となるサブサンプリング位置を求め、求められた前記サブサンプリング位置に最も近似するＮ次元超平面を各パラメータ軸毎に求め、求められた超平面の交点からサブサンプリンググリッド推定位置を求めることが開示されている。
【０００４】
また、これらに関連する技術として、例えば、特許文献２には、画像が撮像装置の特性やノイズなどの影響で変形を受けていても精度が低下しない、サブピクセル精度の位置合わせ手法の提供をすることを課題とし、フィッティング範囲設定手段において、入力画像の微分値の大きな画素の連なりをフィッティング範囲と設定し、関数フィッティング手段において、各フィッティング範囲に対し、単峰性の関数でフィッティングしエッジ位置をサブピクセル精度で求め、次に、ピクセル精度位置合わせ手段において位置合わせを行う。対応候補曲線生成手段において、ある一方の画像上のエッジ点に対して、他方の画像上の対応する複数のエッジ点を曲線で結び対応候補曲線とし、位置ずれ算出手段において、全ての対応候補曲線が交差する点若しくは最も密集する点を求め、その座標を一方の画像からもう一方の画像に対するサブピクセル精度の変位とみなすことが開示されている。
【０００５】
また、例えば、特許文献３には、２枚の画像に基づいてパッチ毎にオプティカルフローを算出した後、オプティカルフローの信頼度に応じて、信頼度の低い領域と信頼度の高い領域とに分け、信頼度の低い領域のオプティカルフローを、その４近傍又は８近傍の領域のオプティカルフローに対してそれらの信頼度に応じた重み付けを行った結果を用いて補間するようにしたオプティカルフロー推定方法が開示されている。
【特許文献１】再表２００４−０６３９９１号公報
【特許文献２】特許第２８９７７７２号公報
【特許文献３】特許第３４３５０８４号公報
【発明の開示】
【発明が解決しようとする課題】
【０００６】
本発明は、このような背景技術の状況の中でなされたもので、演算負荷量を削減して、整数グリッド選択を柔軟にし、かつ、サブピクセル精度の移動位置推定精度を高めるようにした画像処理装置及び画像処理プログラムを提供することを目的としている。
【課題を解決するための手段】
【０００７】
かかる目的を達成するための本発明の要旨とするところは、次の各項の発明に存する。
［１］  複数の画像を入力する画像入力手段と、
前記画像入力手段によって入力された画像を整数グリッド上で変換する変換手段と、
相違度又は類似度を求める関数を近似する２次曲面の係数を算出する係数算出手段と、
前記係数算出手段によって算出された係数のサブ量子化精度での推定を行う推定手段
を具備することを特徴とする画像処理装置。
【０００８】
［２］  前記係数算出手段は、近似の場合に各整数グリッド位置に関する重み付けを行う
ことを特徴とする［１］に記載の画像処理装置。
【０００９】
［３］  前記係数算出手段による重み付けは、各整数グリッドを中心とするボロノイ分割領域の面積又は超体積である
ことを特徴とする［２］に記載の画像処理装置。
【００１０】
［４］  前記ボロノイ分割を行う領域は、２次曲面で近似を行いたい領域内のみに限定する
ことを特徴とする［３］に記載の画像処理装置。
【００１１】
［５］  （５２）式と同等のものが存在するか否かを判断する判断手段
【数５２】
2008209971.tif 000003
をさらに具備し、
前記判断手段によって存在しないと判断された場合には、前記ボロノイ分割を行う領域を広げる
ことを特徴とする［４］に記載の画像処理装置。
【００１２】
［６］  前記領域を広げる場合は、全軸で対称となるように１整数グリッドずつ広げる
ことを特徴とする［５］に記載の画像処理装置。
【００１３】
［７］  整数グリッド上で、整数の範囲で相違度が最小となる点を探索する探索手段
をさらに具備することを特徴とする［６］に記載の画像処理装置。
【００１４】
［８］  前記変換手段による変換は、Ｎ次元であり、Ｎ変数の２次曲面の係数を推定し、
さらに、
（２５）式における行列Ａと同等のものが半正定値であるか否かを検証する検証手段と、
【数２５】
2008209971.tif 000004
Ｍｏｏｒ−Ｐｅｎｒｏｓｅの一般化逆行列と同等のものを算出する一般化逆行列算出手段と、
サブ量子化精度での変換係数を算出する変換係数算出手段
を具備することを特徴とする［１］に記載の画像処理装置。
【００１５】
［９］  前記変換手段による変換は、Ｎ次元であり、Ｎ変数の２次曲面の係数を推定し、
さらに、
（２５）式における行列Ａと同等のものが正定値であるか否かを検証する検証手段と、
【数２５】
2008209971.tif 000005
サブ量子化精度での変換係数を算出する変換係数算出手段
を具備することを特徴とする［１］に記載の画像処理装置。
【００１６】
［１０］  前記変換手段による変換は、平行移動であり、２変数の２次曲面の係数を推定し、
さらに、
（８）式と同等のものが成り立つか否かを検証する検証手段と、
【数８】
2008209971.tif 000006
（９）式と同等のものを算出する算出手段
【数９】
2008209971.tif 000007
を具備することを特徴とする［１］に記載の画像処理装置。
【００１７】
［１１］  さらに、
（１０）式と同等のものを検証する検証手段と、
【数１０】
2008209971.tif 000008
（１２）式と同等のものを算出する算出手段
【数１２】
2008209971.tif 000009
を具備することを特徴とする［８］、［９］又は［１０］に記載の画像処理装置。
【００１８】
［１２］  相違度又は類似度を算出する整数グリッド位置を固定として、予め（５３）式と同等のもの又は（５４）式と同等のものを計算し、（５５）式と同等のもの又は（５６）式と同等のもので２次曲面の係数を推定する
【数５３】
2008209971.tif 000010
【数５４】
2008209971.tif 000011
【数５５】
2008209971.tif 000012
【数５６】
2008209971.tif 000013
ことを特徴とする［１］に記載の画像処理装置。
【００１９】
［１３］  相違度又は類似度を算出する整数グリッド位置は、縦３×横３の９点とし、（３５）式と同等のものと等価の式を用いて、２次曲面の係数を推定する
【数３５】
2008209971.tif 000014
ことを特徴とする［１］に記載の画像処理装置。
【００２０】
［１４］  相違度又は類似度を算出する整数グリッド位置は、縦３×横３の９点とし、（４０）式と同等のものと等価の式を用いて、２次曲面の係数を推定する
【数４０】
2008209971.tif 000015
ことを特徴とする［１］に記載の画像処理装置。
【００２１】
［１５］  さらに、
相違度又は類似度を算出した整数グリッド位置を入力するグリッド位置入力手段と、
行列Ｂを算出する行列算出手段と、
連立一次方程式（１９）式と同等のものを計算する計算手段
【数１９】
2008209971.tif 000016
を具備することを特徴とする［１２］、［１３］又は［１４］に記載の画像処理装置。
【００２２】
［１６］  さらに、
行列Ｂと同等のもののランクを算出する又は（７１）式と同等のものの正則性を検証する検証手段と、
【数７１】
2008209971.tif 000017
相違度を算出する整数グリッドを増加させる整数グリッド増加手段
を具備することを特徴とする［１５］に記載の画像処理装置。
【００２３】
［１７］  さらに、
整数グリッド位置と、行列（５３）式と同等のもの又は行列（５４）式と同等のものの関係を保持する関係保持手段と、
入力された整数グリッド位置が前記関係保持手段に保持されているか否かを判断する判断手段
【数５３】
2008209971.tif 000018
【数５４】
2008209971.tif 000019
を具備し、
前記係数算出手段は、テーブル内の行列（５３）式と同等のもの又は行列（５４）式と同等のものを用いて、２次曲面の係数を算出する
ことを特徴とする［１］に記載の画像処理装置。
【００２４】
［１８］  前記判断手段によって、入力された整数グリッド位置が前記関係保持手段に保持されていないと判断された場合、整数グリッドを増加させる整数グリッド増加手段
をさらに具備することを特徴とする［１７］に記載の画像処理装置。
【００２５】
［１９］  前記推定手段は、画像ブロックの整数移動位置の相違度又は類似度を複数算出し、それを用いてサブピクセル精度の動きベクトルを推定する
ことを特徴とする［１］に記載の画像処理装置。
【００２６】
［２０］  さらに、
全ての変換係数を０．５グリッド分ずらす画像を作成する画像作成手段
を具備し、
前記推定手段は、前記画像作成手段によって生成された画像で求めた位置と、ずらす前の画像で求めた位置とから、（４５）式と同等のもので位置を求める
【数４５】
2008209971.tif 000020
ことを特徴とする［１］に記載の画像処理装置。
【００２７】
［２１］  さらに、
一つ以上の変換係数を０．５グリッド分ずらす画像を作成する画像作成手段
を具備し、
前記推定手段は、前記画像作成手段によって生成された画像で求めた位置と、ずらす前の画像で求めた位置とから、（４７）式と同等のもの、（４８）式と同等のもの又は（４９）式と同等のものを用いた式で位置を求める
【数４７】
2008209971.tif 000021
【数４８】
2008209971.tif 000022
【数４９】
2008209971.tif 000023
ことを特徴とする［１］に記載の画像処理装置。
【００２８】
［２２］  画像処理システムを、
複数の画像を入力する画像入力手段と、
前記画像入力手段によって入力された画像を整数グリッド上で変換する変換手段と、
相違度又は類似度を求める関数を近似する２次曲面の係数を算出する係数算出手段と、
前記係数算出手段によって算出された係数のサブ量子化精度での推定を行う推定手段
として機能させることを特徴とする画像処理プログラム。
【発明の効果】
【００２９】
本発明にかかる画像処理装置及び画像処理プログラムによれば、本構成を有していない場合に比較して、演算負荷量を削減して、整数グリッド選択を柔軟にし、かつ、サブピクセル精度の移動位置推定精度を高めることができる。
【発明を実施するための最良の形態】
【００３０】
＜１．前提＞
本実施の形態は、２枚以上の画像の位置合わせをサブピクセル精度（あるいはサブ量子化インデクス精度）で行うものである。
まず、本実施の形態を説明する前に、特許文献１に記載されている技術を中心にして画像の位置合わせに関する技術の説明を行う。
特許文献１に記載されている技術の前提を説明する。以下に示す前提は本実施の形態の前提でもある。
【００３１】
２枚の画像ｆ（ｕ，ｖ）と、画像ｇ（ｕ，ｖ）が入力されたとする。また、画像ｆ（ｕ，ｖ）内の領域Ｗを定義する。
画像ｇ（ｕ，ｖ）をｕ方向にｘ、ｖ方向にｙ移動させると、ｇ（ｕ−ｘ，ｖ−ｙ）となる。このとき、領域Ｗ内における画像ｆ（ｕ，ｖ）と画像ｇ（ｕ−ｘ，ｖ−ｙ）の差（相違度）ができるだけ小さくなるようにしたい。このような移動量（ｘ，ｙ）を求めることが特許文献１に記載の技術（画像位置合わせ技術）の目的である。
【００３２】
相違度Ｒは移動量（ｘ，ｙ）の関数Ｒ（ｘ，ｙ）として表すことができる。
相違度関数の例としては、画素値の２乗誤差、絶対誤差が考えられる。
画素値の２乗誤差を用いる場合の相違度関数ＲＳＳＤ（ｘ，ｙ）は（１）式のように表すことができる。ただし、ＳＳＤは、Ｓｕｍ  ｏｆ  Ｓｑｕａｒｅｄ  Ｄｉｆｆｅｒｅｎｃｅｓの略である。
【数１】
2008209971.tif 000024
【００３３】
画素値の絶対誤差を用いる場合の相違度関数ＲＡＳＤ（ｘ，ｙ）は（２）式のように表すことができる。ただし、ＡＳＤは、Ｓｕｍ  ｏｆ  Ａｂｓｏｌｕｔｅ  Ｄｉｆｆｅｒｅｎｃｅｓの略である。
【数２】
2008209971.tif 000025
あるいは差（相違度）を小さくするということは、類似度を大きくするということであるとみなすこともできる。この場合、類似度としては、正規化相関係数が考えられる。相違度ではなくて類似度を用いる場合は、類似度を最大化する移動量を求めればよい。以下の説明では相違度を例として説明し、相違度を最小化する方法を述べるが、類似度を用いる場合には、単に最小化を最大化と読み替えればよい。
【００３４】
画素値の正規化相関係数を用いる場合の相違度関数ＲＺＮＣＣ（ｘ，ｙ）は（３）式のように表すことができる。ただし、ＺＮＣＣは、Ｚｅｒｏ−ｍｅａｎ  Ｎｏｒｍａｌｉｚｅｄ  Ｃｒｏｓｓ−Ｃｏｒｒｅｌａｔｉｏｎの略である。
【数３】
2008209971.tif 000026
【００３５】
これらのどの相違度（あるいは類似度）を用いてもよい。そこで、以下、ＲＳＳＤ（ｘ，ｙ）を最小化する（ｘ，ｙ）の求め方を示す。以下、ＲＳＳＤ（ｘ，ｙ）を単にＲ（ｘ，ｙ）と記述する。
特許文献１に記載の技術は、ｘ，ｙを整数以下の精度まで求める手法を示している。
特許文献１に記載の技術の前提として、ｘ，ｙが整数の精度では最適な位置を求めることができているとする。
このような整数精度で最適化されたｘ，ｙを、Ｘ０，Ｙ０とする。
特許文献１には明記されていないが、このようなｘ，ｙの組み合わせを発見することは可能である。一例を挙げると、あらゆる整数の（ｘ，ｙ）の組み合わせに関してＲ（ｘ，ｙ）を調べて、Ｒが最小となる（ｘ，ｙ）を採用すればよい。（ｘ，ｙ）の範囲は、画像ｆと画像ｇが重なりを持つ範囲内に限定されるため、有限である。
あるいは、階層的に行うなどの種々の高速化手法を用いてもよい。
【００３６】
以上の前提に立つことで、課題を下記のように書き換えることができる。
課題：「Ｒ（Ｘ０＋ｘ，Ｙ０＋ｙ）を最小化する。ただし、−１＜ｘ＜１、かつ、−１＜ｙ＜１」
ここで、Ｘ０，Ｙ０は、単なるオフセットである。毎回記述することは煩雑となるため、ここでは画像ｆと画像ｇは既に整数精度で移動位置が最適化されているものとみなす。このようにすると、Ｘ０＝０，Ｙ０＝０とできるため、前記課題を、
課題：「Ｒ（ｘ，ｙ）を最小化する。ただし、−１＜ｘ＜１、かつ、−１＜ｙ＜１」
のように簡易に記述することが可能である。以下、この前提を用いる。
【００３７】
＜２．特許文献１に記載の技術＞
特許文献１に記載の技術の概略は、以下の通りである。
（１）ｘ，ｙが整数の場合のＲの値をまず求める。以下、ｘ，ｙが整数となるような位置を整数グリッドと称する。
（２）上記整数グリッド位置のＲの値を用いて、Ｒが最小となるｘ，ｙ（−１＜ｘ＜１、−１＜ｙ＜１）を求める。
ここで、用いる整数グリッドを図１６に示す。
特許文献１に記載の技術では、まずｘ＝０の直線において、Ｒを最小とするｙの値を推定する。これをｙ（０）とする。特許文献１に記載の技術では、Ｒ（０，−１）、Ｒ（０，０）、Ｒ（０，１）の３点を通る放物線関数を決定して、その放物線関数を最小とする位置をｙ（０）とする。
次に、ｓ＝１の直線において、Ｒを最小とするｙの値を推定する。これをｙ（１）とする。ｙ（１）を求めるとき、まず、Ｒ（１、−３）〜Ｒ（１，３）の中で最小となる点を発見する。その点を中心とした３点のＲの値を抽出して、その３点を通る放物線関数を決定する。さらに、その放物線関数を最小とする位置をｙ（１）とする。
次にｘ＝−１の直線において、Ｒを最小とするｙの値を推定する。これをｙ（−１）とする。ｙ（−１）の求め方はｙ（１）と同様である。
【００３８】
さらに、ｙ＝０の直線においてＲを最小とするｘの値ｘ（０）、ｙ＝１の直線においてＲを最小とするｘの値ｓ（１）、及び、ｙ＝１の直線においてＲを最小とするｘの値ｘ（−１）を求める。
さらに、
直線１  （０，ｙ（０）），（１，ｙ（１）），（−１，ｙ（−１））を通る直線
直線２  （ｘ（０），０），（ｘ（１），１），（ｘ（−１）．−１）を通る直線
を求める。上記３点が直線上に並んでいない場合、最小２乗法で直線を決める。
さらに、直線１と直線２の交点を求める。この交点が求めるサブピクセル精度位置（ｘ，ｙ）である。
特許文献１に記載の技術を用いて推定したサブピクセル精度移動位置を図１７に示す。つまり、図１７に示した直線１と直線２とが交わる点が特許文献１に記載の技術が推定した移動位置ということになる。
【００３９】
＜３．実施の形態が解決しようとする課題＞
特許文献１に記載の技術では、以下に示すような課題が存在する。
（１）ｘ，ｙ＝±３の範囲の計２５整数グリッドのＲの値を取得しているため、演算処理量が多い。
つまり、１グリッドのＲの値を算出するためには、（１）式に示されるように、領域Ｗの画素数分だけ演算を行わなければならないため、演算負荷が大きい。
（２）ｘ＝±１又はｙ＝±１の場合に、最小値を求める負荷がかかる。
（３）ｘ，ｙ＝±３の範囲に限定しているため、最小値がこの範囲外にある場合、推定精度が低下する。
（４）各直線上（ｘ＝０，ｘ＝１，ｘ＝−１，ｙ＝０，ｙ＝−１，ｙ＝１等）で、最小値を求めるためには、各直線上で最低３点のＲの値が必要である（放物線を推定するには３点が必要である）。つまり、整数グリッドの選択基準に柔軟性がない。
（５）相違度空間Ｒが２次元ガウス関数で近似できる場合に前述のような直線の交点が最適であることが示されているが、実際の画像でこのような近似が成り立つ保証はなく、精度が落ちる虞がある。
そこで、本実施の形態は、特許文献１に記載の技術よりも演算負荷量を削減して、かつ、整数グリッド選択を柔軟にし、かつ、サブピクセル精度の移動位置推定精度を高めるようにするものである。
【００４０】
＜４．実施の形態＞
＜４．１．実施の形態の基本的アーキテクチャ＞
＜４．１．１．２次元移動位置＞
本実施の形態は、相違度空間Ｒにおける相違度を２変数の２次曲面にマッピングし、その２次曲面を最小化する移動量を求めるものである。
すなわち、いま、（５７）式として、相違度Ｒをｘとｙの（４）式に示す関数（２次曲面）で近似する。
【数５７】
2008209971.tif 000027
【数４】
2008209971.tif 000028
【００４１】
この関数Ｒを最小化するｘ，ｙの必要条件は、δＲ／δｘ＝０、及び、δＲ／δｙ＝０であるから、（５）式のようになる。
【数５】
2008209971.tif 000029
結局、（６）式を解けばよいことになる。
【数６】
2008209971.tif 000030
この（６）式を以下の条件１〜４に分けて解く。
【００４２】
（条件１）
２次曲面（４）が最小値を１点だけ持つためには、行列式である（７）式が正定値（ｐｏｓｉｔｉｖｅ  ｄｅｆｉｎｉｔｅ）であればよい。
【数７】
2008209971.tif 000031
すなわち、（８）式が成り立てばよい。
【数８】
2008209971.tif 000032
このときは、単純に（６）式を解けばよい。すなわち、（９）式のようになる。
【数９】
2008209971.tif 000033
なお、量子化精度は、１画素であってもよく、２画素以上であってもよい。
【００４３】
（条件２）
次に、Ａの行列式が０で、かつ最小値が存在する場合を考える。この場合は、（１０）式となる。
【数１０】
2008209971.tif 000034
このとき、行列Ａが正則ではないので、Ｓの最小値を与えるｔが無限に存在することになる（（６）式が不定の方程式となる）。ここでは、ｔのノルムが最小のものを採用する。なぜならば、整数に限定した場合に（０，０）が最小となるようにしているため、これは妥当な考え方である。つまり、直線の（１１）式の中で、ノルムが最小の（ｘ，ｙ）を求める。
【数１１】
2008209971.tif 000035
解は、（１２）式のようになる。
【数１２】
2008209971.tif 000036
【００４４】
（条件３）
次に、ａ＝ｂ＝ｃ＝０、すなわち、Ａ＝Ｏ（ゼロ行列）のときは、（５８）式とすればよい。
【数５８】
2008209971.tif 000037
【００４５】
（条件４）
その他の場合、最小値を持つことはない。この場合、（５８）式を採用する。
【００４６】
ここで、Ｍｏｏｒｅ−Ｐｅｎｒｏｓｅの一般化逆行列を用いて前記（６）式を解く方法を考慮する。係数行列のＭｏｏｒｅ−Ｐｅｎｒｏｓｅの一般化逆行列を用いて線形方程式を解く場合、かつ、係数行列が正則ではない場合、ノルムが最小の解を求めることになる。よって、Ｍｏｏｒｅ−Ｐｅｎｒｏｓｅの一般化逆行列を用いることによって、上記の（条件１）〜（条件３）を包含することが可能となる。
【００４７】
以上をまとめる。
ａ≧０、かつ、ｂ≧０、かつ、ａｂ−ｃ２≧０のとき（（条件１）〜（条件３）の場合）、（１３）式として前記（６）式を解けばよい。
【数１３】
2008209971.tif 000038
それ以外のとき（条件４の場合）、（５８）式とする。
【００４８】
さらに、計算で得られた（５９）式が、ｘ＞１、ｘ＜−１、又は、ｙ＞１、ｙ＜−１となる場合がありうる。このときは、安定化のため（５８）式としてもよいと考えられる。
【数５９】
2008209971.tif 000039
ただし、前述した説明では、前記（４）式の係数ａ〜ｅを示していない。
【００４９】
次の問題は、前記（４）式の係数ａ〜ｅの推定である（ｆは前記（６）式の中にないため、不要）が、これは以下のようにする。
まず、（６０）式を整数グリッドであるとする。
【数６０】
2008209971.tif 000040
これら複数個の整数グリッドを用いて、相違度Ｒ（ｘｉ，ｙｉ）を計算する。以下、表示を簡単にするため、Ｒｉ＝Ｒ（ｘｉ，ｙｉ）とする。
以下、（１４）式とする。
【数１４】
2008209971.tif 000041
このとき、Ｒｉは、（１５）式のように書ける。
【数１５】
2008209971.tif 000042
そこで、（１６）式とする。複数の点（ｉ＝０，１，２．．．）で（１５）式を求めることができるので、（１７）式とできる。
【数１６】
2008209971.tif 000043
【数１７】
2008209971.tif 000044
ここで、（１８）式とすれば、連立１次方程式（１９）式を形成することができる。
【数１８】
2008209971.tif 000045
【数１９】
2008209971.tif 000046
ａは未知数６個のベクトルであるため、行列Ｂのランクが６以上であれば解くことができる。整数グリッドの位置は、行列Ｂのランクが６以上であれば何でもよい。このように、本実施の形態では、整数グリッドの位置を非常に自由に取ることができる。
一般的には、整数グリッドが２次元状に分布するように、整数グリッドの個数を６以上取ればよい。このとき行列Ｂの行数は６以上となるので、最小２乗法で解くことができる。
【数２０】
2008209971.tif 000047
（２０）式において、ａの最終要素ｆを算出する必要はないため、最終行の乗算は行わなくともよい。
【００５０】
＜４．１．２．重み付け＞
さらに、（２０）式では、各画素位置ｉに関する重みを一定として最小２乗法を適用したが、各整数グリッド位置ｉに関する重み付けを行うことも可能である。
整数グリッド位置ｉに関する重みをｗｉ＞０とする。
さらに、（２１）式とすると、最小２乗近似は、（２２）式とすればよい。
【数２１】
2008209971.tif 000048
【数２２】
2008209971.tif 000049
行列Ｗが単位行列である場合、すなわちＷ＝Ｉ（単位行列）の場合、重み付けを行わないときに一致する。
以下、重み係数設定の一つの手法を示す。
まず、各整数グリッドの値は相違度空間を代表していると考える。代表している相違度空間の面積（又は超体積）を重み係数とすればよいと考える。
ここでは、各整数グリッドを中心として相違度空間を領域πの中で、ボロノイ分割する。分割された各領域の面積（又は超体積）を重み係数とする。
ボロノイ分割とは、近い整数グリッド間の垂直２等分線で空間分割を行う方法である。Ｎ次元空間であれば、垂直超２等分平面で分割することになる。
【００５１】
ここで、相違度空間の範囲は、２次曲面で近似を行いたい範囲に限定することで精度を高めることができる。例えば、小数点以下の画素位置を推定したい場合、２次曲面で近似を行いたい範囲を−１≦ｘ≦１、−１≦ｙ≦１とする。この範囲をπとする。πを、−１≦ｘ≦１、−１≦ｙ≦１の領域を考える。この領域内において、図１８に示すように９点（図１８に示す格子上の０〜８が該当）の整数グリッドがあるとする。この９点を中心とするボロノイ分割結果は、図１９に示すようになる。
図１９から分かるように、
（１）整数グリッド０を中心とする領域の面積は、１
（２）整数グリッド１、２、３、４を中心とする領域の面積は、１／２
（３）整数グリッド５、６、７、８を中心とする領域の面積は、１／４
となる。この結果より、図２０に示すように重み係数を設定するとすると、ｗａ＝０．５，ｗｂ＝０．２５となる。
【００５２】
このようにボロノイ分割して重み係数を決定すると、重み付けを行わない場合（すなわちＷ＝Ｉの場合）には計算ができる場合であっても、重みが０となってしまって、（６１）式が正則にならない場合がありうる。例えば、πが−１≦ｘ≦１、−１≦ｙ≦１の領域で、整数グリッドが、（０，０）、（０，２）、（０，−２）、（２，０）、（−２，０）、（２，２）、（−２，２）、（２，−２）、（−２，−２）の場合（図２１参照）中心以外の重みが０となるため、上記のようなことがありうる。
【数６１】
2008209971.tif 000050
このときは領域πを−２≦ｘ≦２、−２≦ｙ≦２に広げてやればよい。
例えば、計画性なくπの領域を広げると、精度が悪くなるため、（６１）式が正則であるという条件下で、できるだけπは小さな領域に限定することが望ましい。そのためには、領域πを−ｐ≦ｘ≦ｐ、−ｐ≦ｙ≦ｐとして、ｐを１から順に、２、３と増加させていき、（６１）式が正則になった時点で終了すればよい。
なお、重み係数設定方法はこれに限らない。予め固定しておいてもよい。あるいは、整数グリッド位置と重み係数の関係をテーブルに保持しておいてもよい。あるいは、整数グリッド位置と前記（５４）式の関係をテーブルに保持しておいてもよい。もちろん、テーブル構造に限る必要はなく、リンク構造等であってもよく、これらの対応関係を記憶できる記憶手段であればよい。以下、記憶手段としてテーブル構造を主として説明する。
【００５３】
＜４．１．３．多次元への拡張＞
前述の説明では、ｆとｇの関係を平行移動に限った場合であり、求めるベクトルは２次元であった。ここでは、これをＮ次元へ拡張する。
ｆとｇの間の変換は、平行移動以外にも様々なものが考えられる。その際に、変換を表現する変換係数をベクトルとしてみなせば、前節までの平行移動ベクトルと同様に考えることができる。
例えば、回転を含めれば３次元となる。相似変換（回転＋相似拡大）では４次元、縦横の拡縮を入れると５次元、アフィン変換では６次元、射影変換では８次元のベクトルを求めればよい。
２次元の場合に求めたベクトルｔをここでは、Ｎ次元ベクトルｈに拡張する。ｈは、射影変換行列の各要素をベクトル化したものと捉えてよい。
回転と移動の３次元の変換であれば、ｈは、２次元の移動ベクトルと、回転角を要素として持つと考えればよい。ｈの各要素は、擬似的に量子化された整数グリッドを持つものとする。例えば回転角であれば、「１０度を１グリッドとして量子化する」などと定めればよい。
ちなみに、同様に平行移動に関しても、１画素を１グリッドとする必然性はないので、任意のｐ画素を１グリッドとすることができる。「任意のｐ画素を１グリッドとすることができる。」というのは、前節までの２次元の場合でも成り立つ。
【００５４】
一般には、ベクトルｈは、射影変換行列の各要素をベクトル化したものと捉えてよい。いま、（６２）式とする。
【数６２】
2008209971.tif 000051
まず、２次元の場合の例を前記ベクトルｈを用いて記述しなおす。（６３）式を用いて、（２３）式のように変形することができる。
【数６３】
2008209971.tif 000052
【数２３】
2008209971.tif 000053
２次元のときは、（２３）式のように書くことができる。これをそのままＮ次元に拡張する。すなわち、Ｒを（２４）式の２次曲面で近似する。
【数２４】
2008209971.tif 000054
ただし、行列ＡはＮ×Ｎの対称行列、ｄはＮ次元のベクトル、ｆはスカラー値とできる。行列Ａ、ベクトルｄ、ベクトルｈは（２５）式で与えられる。
【数２５】
2008209971.tif 000055
Ｒを最小化する必要条件は、（２６）式とすればよい。
【数２６】
2008209971.tif 000056
Ａは対称行列なので、（２７）式とできる。
【数２７】
2008209971.tif 000057
【００５５】
この（２７）式を解くには、
（１）行列Ａが半正定値（Ａの固有値が全て０以上）のとき、（２８）式とすればよい。
【数２８】
2008209971.tif 000058
ただし、（２８）式内のＡ＋は、Ｍｏｏｒｅ−Ｐｅｎｒｏｓｅの一般化逆行列である。
（２）前記（１）以外のとき、（２９）式とすればよい。
【数２９】
2008209971.tif 000059
さらに、
（３）あるｉで、ｈｉ＞１  ｏｒ  ｈｉ＜−１のとき、ｈｉ＝０としてもよい。
【００５６】
次に、行列Ａとベクトルｄの算出方法を述べる。
未知ベクトルａを（３０）式のように定義する。
【数３０】
2008209971.tif 000060
（３０）式は、Ｎ（Ｎ＋１）／２＋Ｎ＝Ｎ（Ｎ＋３）／２次元ベクトルとなる。
前記（２４）式を（３０）式のベクトルａを用いて記述しなおすと、（３１）式となる。
【数３１】
2008209971.tif 000061
【００５７】
ここで、整数グリッド番号をｉとして、番号ｉの整数グリッドのときの相違度をＲｉとする。また、整数グリッドｉのときのベクトルｈを、（６４）式とする。
【数６４】
2008209971.tif 000062
各整数グリッド点毎に（３１）式を得ることができるため、（３２）式として、前記（１８）式を用いて、Ｂとｒを求めて、前記（２０）式又は前記（２２）式を適用すれば、ベクトルａを求めることができる。
【数３２】
2008209971.tif 000063
【００５８】
＜４．２．具体的演算例＞
＜４．２．１．重み付けなし＞
前述した説明に関して、具体的な数値例を示す。
ここでは２次元の例を示す。
また、相違度を算出する整数グリッドを（６５）式に示す９つとする。
【数６５】
2008209971.tif 000064
図２２は、（６５）式に示す整数グリッドを図示したものである。
この整数グリッドの場合に、Ｂは（３３）式となる。
【数３３】
2008209971.tif 000065
また、（６６）式となる。
【数６６】
2008209971.tif 000066
ここで、前記（５３）式を計算すると、（３４）式のようになる。
【数３４】
2008209971.tif 000067
よって、係数ベクトルａは、（３５）式で求めればよい。
【数３５】
2008209971.tif 000068
【００５９】
ｆは算出の必要がないため、（３５）式の最終行の乗算は省くことができる。（３６）式に示すように、式が簡単になる。
【数３６】
2008209971.tif 000069
この（３６）式を利用して、以下の手順でｔを求める。
（１）ａ＞０、ｂ＞０、ａｂ−ｃ２＞０のとき、（３７）式を用いる。
【数３７】
2008209971.tif 000070
（２）ａ＞０、ｂ＞０、ａｂ−ｃ２＝０のとき、前記（１２）式を用いる。
（３）上記以外のとき、ｔ＝（０，０）Ｔとする。
【００６０】
＜４．２．２．重み付け有り＞
次に重み付けを行う場合の例を示す。グリッド点は前節と同じとする。すなわち行列Ｂは同じである。また、前記（６６）式となる。
ここでは、重み付け行列Ｗを（３８）式のように決める。
【数３８】
2008209971.tif 000071
例えば、ｗａ＝０．５，ｗｂ＝０．２５とすると、（３９）式のように計算できる。
【数３９】
2008209971.tif 000072
前節と同様にｆが算出不要なので、（４０）式のようになる。
【数４０】
2008209971.tif 000073
（４０）式に示されるように、非常に単純な演算となる。このａ〜ｅを用いて前節と同様にｔを計算すればよい。
さらに、（４０）式では、行列の要素に０が存在しているため、全てを計算する必要はない。よって、（４１）式、（４２）式、（４３）式、（４４）式に示すように書きなおすこともできる。
【数４１】
2008209971.tif 000074
【数４２】
2008209971.tif 000075
【数４３】
2008209971.tif 000076
【数４４】
2008209971.tif 000077
【００６１】
ここで、この手法を用いて、２次曲面を推定した例を示しておく。
例えば、一般的な実験用画像である画像Ｂａｒｂａｒａの１／１０縮小画像を用いて人工的に（−０．３，−０．２）画素ずれた画像を作成した。この画像をずらして、真のＳＳＤ値を測定した。
この場合、ＳＳＤの値は、図２３に示すように変化する。
さらに、前記画像を１画素単位でずらして求めた９つの整数グリッド点情報を用いて、２次曲面を推定した。
近似２次曲面は図２４に示すようになる。推定結果を破線で示す。
図２４に示されるように、真のＳＳＤ値と推定２次曲面がほとんど一致していることが分かる。
この場合の最小を与える（ｘ，ｙ）の推定位置は、（−０．２７８，−０．１８６）となる。真の値は（−０．３，−０．２）であり、実用上十分な精度が得られていることが分かる。
さらに、分かりやすいように等ＳＳＤ線を図２５に示す。
【００６２】
＜４．３．高精度化方式＞
特許文献１には、放物線でＲを近似したときに発生する位置誤差を解消する方式の例が示されている。
これは、放物線による近似では、画素の移動距離が整数グリッドと±０．２５の差があるときに最も誤差が大きくなる現象を回避するものである。具体的には、０．５画素ずらした画像を（例えば線形補間で）作成し、その画像を用いて推定を行う。この推定誤差は、ずらさない画像の推定誤差と逆位相となるため、双方の結果の平均値を取ることによって、誤差をキャンセルする。
この高精度化方式を適用する場合、従来例では、１次元の推定をまず行うため、ｘ方向に０．５画素ずれた画像を用いてｘ方向の誤差キャンセルを行う。また、ｙ方向に０．５画素ずれた画像を用いてｙ方向の誤差キャンセルを行う。すなわち、Ｎ次元パラメタを推定する場合、Ｎ個の０．５画素ずれた画像を作成する必要がある。
本実施の形態に関しても、この高精度化方式が適用可能である。ただし、本手法ではＮ次元を一気に推定するため、１個だけ（全てのパラメタに関して０．５インデクスずれた）画像を作成して、その画像を用いて推定した結果との平均値を取ればよい。特許文献１に記載の技術と比較して、補間画像作成数を１／Ｎとすることが可能である。この手法を「高精度化１手法」とする。
全てのｈｉ（ｉ＝１，２，．．．，Ｎ）に関して０．５グリッド（画素）ずれた画像をｇから線形補間等で作成する。この画像を用いて位置合わせを行った結果を、ｈｓとする。また、元々のｇ画像を用いて求めたｈをｈｏとする。最終的に求めるｈは、（４５）式のようになる。もちろん、ｈｓは、０．５グリッドのシフト分を補償済みの値とする。
【数４５】
2008209971.tif 000078
【００６３】
あるいは、本実施の形態でも１次元のＮ方向の補間画像を用いた高精度化も考えられる。複数の補間結果を用いることができるため、さらに精度を上げることができる。
あるｉに対してのみ０．５グリッド（画素）ずれた画像をｇから線形補間等で作成する。この画像を用いて位置合わせを行った結果を、（６７）式とする。
【数６７】
2008209971.tif 000079
また、元々のｇ画像を用いて求めたｈを（６８）式とする。
【数６８】
2008209971.tif 000080
このとき、最終的に求める（６９）式は、（４６）式として求められる。
【数６９】
2008209971.tif 000081
【数４６】
2008209971.tif 000082
これは、「シフトしていない画像の要素の平均値」と、「シフトした画像の要素」との平均値を求めるものである。２次元の場合は、（４７）式となる。この手法を「高精度化２手法」とする。
【数４７】
2008209971.tif 000083
【００６４】
高精度化２手法では、ベクトルの各要素毎に、「シフトしていない画像の平均値」と、「シフトした画像の平均値」を求めていたが、平均値を求めるときには、「シフトした画像」と、「その他の画像の要素を適当に選択したとき」との平均値であれば何でもよい。つまり、（４８）式、（４９）式、（５０）式のようになる。
（４８）式、（４９）式で関数ａｖｅｒａｇｅｊ（）は、変数ｊを変化させて平均値を計算した結果を示す。つまり、ｉは変化させない。また、ｊの変化の範囲は様々にすることができる。
【数４８】
2008209971.tif 000084
【数４９】
2008209971.tif 000085
【数５０】
2008209971.tif 000086
【００６５】
以下、図面に基づき本発明を実現するにあたってのより具体的な各種の実施の形態を説明する。
＜４．４．第１の実施の形態＞
図１は、第１の実施の形態の概念的なモジュール構成図を示している。
なお、モジュールとは、一般的に論理的に分離可能なソフトウェア、ハードウェア等の部品を指す。したがって、本実施の形態におけるモジュールはプログラムにおけるモジュールのことだけでなく、ハードウェア構成におけるモジュールも指す。それゆえ、本実施の形態は、プログラム、システム及び方法の説明をも兼ねている。また、モジュールは機能にほぼ一対一に対応しているが、実装においては、１モジュールを１プログラムで構成してもよいし、複数モジュールを１プログラムで構成してもよく、逆に１モジュールを複数プログラムで構成してもよい。また、複数モジュールは１コンピュータによって実行されてもよいし、分散又は並列環境におけるコンピュータによって１モジュールが複数コンピュータで実行されてもよい。なお、１つのモジュールに他のモジュールが含まれていてもよい。また、以下、「接続」とは物理的な接続の他、論理的な接続（データの授受、指示等）を含む。
また、システムとは、複数のコンピュータ、ハードウェア、装置等がネットワーク等の通信手段で接続されて構成されるほか、１つのコンピュータ、ハードウェア、装置等によって実現される場合も含まれる。
【００６６】
まずは、２次元で整数グリッドが固定の場合の例を示す。
第１の実施の形態では、図２２に示された９点の固定の整数グリッドを用いる。既に整数の範囲では中心の整数グリッドにおける相違度が最小であるとする。
この例では、重みを付けていてもよいし、重みを付けていなくてもよい。
【００６７】
第１の実施の形態は、図１に示すように、整数グリッド上の相違度算出モジュール１１、２次曲面係数算出モジュール１２、位置算出モジュール１３を有している。整数グリッド上の相違度算出モジュール１１は２次曲面係数算出モジュール１２に、２次曲面係数算出モジュール１２は整数グリッド上の相違度算出モジュール１１、位置算出モジュール１３に、位置算出モジュール１３は２次曲面係数算出モジュール１２に、接続されている。
まず、整数グリッド上の相違度算出モジュール１１は、入力画像１及び入力画像２を入力する。画像を入力するとは、具体的には、スキャナによって画像を入力すること、ファックスによって通信回線を介して外部機器から画像を入力すること、記憶媒体から画像を入力すること等がある。
整数グリッド上の相違度算出モジュール１１では、画像１の領域Ｗを整数範囲で動かして、整数グリッド上の９点の相違度を算出する。ここで、既に整数の範囲では中心の整数グリッドにおける相違度が最小であるとする。算出方法は例えば前記（１）式を用いる。出力は前記（１８）式に示されるベクトルｒである。
２次曲面係数算出モジュール１２では、２次曲面係数ａ，ｂ，ｃ，ｄ，ｅを算出する。例えば、前記（４０）式を用いればよい。重みを付けない場合は、前記（３６）式を利用すればよい。
位置算出モジュール１３では、まず前記（８）式が成り立つかどうかを検証する。
前記（８）式が成り立つ場合、前記（９）式を用いて、最終的な位置ｔを求める。
前記（８）式が成り立たない場合、かつ、前記（１０）式が成り立つ場合、前記（１２）式を用いて、最終的な位置ｔを求める。
前記以外の場合、（７０）式とする。
【数７０】
2008209971.tif 000087
あるいは、以下のように簡単化してもよい。
前記（８）式が成り立つ場合、前記（９）式を用いて、最終的な位置ｔを求める。
前記以外の場合、（７０）式とする。
【００６８】
＜４．５．第２の実施の形態＞
第２の実施の形態は、２次元で整数グリッドが可変の場合に対応できるものである。
第１の実施の形態では、前記（５３）式を予め計算しておくことができたが、整数グリッド位置が可変の場合は、予め計算しておくことができない。そこで、毎回、行列Ｂを算出するようにしたものである。
第２−１の実施の形態は、図２に示すように、整数グリッド上の相違度算出モジュール２１、２次曲面係数算出モジュール２２、位置算出モジュール２３、行列Ｂ算出モジュール２４、（５３）式算出モジュール２５を有している。整数グリッド上の相違度算出モジュール２１は２次曲面係数算出モジュール２２、行列Ｂ算出モジュール２４に、２次曲面係数算出モジュール２２は整数グリッド上の相違度算出モジュール２１、（５３）式算出モジュール２５、位置算出モジュール２３に、位置算出モジュール２３は２次曲面係数算出モジュール２２に、行列Ｂ算出モジュール２４は整数グリッド上の相違度算出モジュール２１、（５３）式算出モジュール２５に、（５３）式算出モジュール２５は２次曲面係数算出モジュール２２、行列Ｂ算出モジュール２４に、それぞれ接続されている。
整数グリッド上の相違度算出モジュール２１は、処理結果である整数グリッド位置を行列Ｂ算出モジュール２４へ渡す以外は、第１の実施の形態の整数グリッド上の相違度算出モジュール１１と同様である。
２次曲面係数算出モジュール２２は、（５３）式算出モジュール２５からの処理結果を受け取る以外は、第１の実施の形態の２次曲面係数算出モジュール１２と同様である。
位置算出モジュール２３は、第１の実施の形態の位置算出モジュール１３と同様である。
行列Ｂ算出モジュール２４は、行列Ｂを、前記（１６）式、前記（１８）式を用いて算出する。
（５３）式算出モジュール２５は、通常の行列演算手法を用いて、前記（５３）式を算出する。
【００６９】
又は、図３に示すように、整数グリッド位置から重み行列Ｗを算出し、さらに前記（５４）式を算出してもよい。
図３に示す第２−２の実施の形態は、整数グリッド上の相違度算出モジュール３１、２次曲面係数算出モジュール３２、位置算出モジュール３３、行列Ｂ算出モジュール３４、重み行列Ｗ算出モジュール３５、（５４）式算出モジュール３６を有している。整数グリッド上の相違度算出モジュール３１は２次曲面係数算出モジュール３２、行列Ｂ算出モジュール３４、重み行列Ｗ算出モジュール３５に、２次曲面係数算出モジュール３２は整数グリッド上の相違度算出モジュール３１、（５４）式算出モジュール３６、位置算出モジュール３３に、位置算出モジュール３３は２次曲面係数算出モジュール３２に、行列Ｂ算出モジュール３４は整数グリッド上の相違度算出モジュール３１、（５４）式算出モジュール３６に、重み行列Ｗ算出モジュール３５は整数グリッド上の相違度算出モジュール３１、（５４）式算出モジュール３６に、（５４）式算出モジュール３６は行列Ｂ算出モジュール３４、重み行列Ｗ算出モジュール３５、２次曲面係数算出モジュール３２に、それぞれ接続されている。
整数グリッド上の相違度算出モジュール３１は、重み行列Ｗ算出モジュール３５にも処理結果である整数グリッド位置を渡す以外は、第２−１の実施の形態の整数グリッド上の相違度算出モジュール２１と同様である。
２次曲面係数算出モジュール３２、位置算出モジュール３３、行列Ｂ算出モジュール３４は、それぞれ第２−１の実施の形態の２次曲面係数算出モジュール２２、位置算出モジュール２３、行列Ｂ算出モジュール２４と同様である。
（５４）式算出モジュール３６は、重み行列Ｗ算出モジュール３５から処理結果である重み行列Ｗを受け取る以外は、第２−１の実施の形態の（５３）式算出モジュール２５と同様である。
重み行列Ｗ算出モジュール３５は、整数グリッド位置から重み行列Ｗを算出し、その重み行列Ｗを（５４）式算出モジュール３６へ渡す。
【００７０】
＜４．６．第３の実施の形態＞
前述の第２の実施の形態は、前記（７１）式が正則であることが保証されている場合である。前記（７１）式が正則であることが保証されていない場合は第３の実施の形態のような構成をとる必要がある。
第３の実施の形態は、図４に示すように、整数グリッド上の相違度算出モジュール４１、２次曲面係数算出モジュール４２、位置算出モジュール４３、行列Ｂ算出モジュール４４、（７１）式正則性検証モジュール４５、（５３）式算出モジュール４６を有している。整数グリッド上の相違度算出モジュール４１は２次曲面係数算出モジュール４２、行列Ｂ算出モジュール４４、（７１）式正則性検証モジュール４５に、２次曲面係数算出モジュール４２は整数グリッド上の相違度算出モジュール４１、（５３）式算出モジュール４６、位置算出モジュール４３に、位置算出モジュール４３は２次曲面係数算出モジュール４２に、行列Ｂ算出モジュール４４は整数グリッド上の相違度算出モジュール４１、（７１）式正則性検証モジュール４５に、（７１）式正則性検証モジュール４５は整数グリッド上の相違度算出モジュール４１、行列Ｂ算出モジュール４４、（５３）式算出モジュール４６に、（５３）式算出モジュール４６は（７１）式正則性検証モジュール４５、２次曲面係数算出モジュール４２に、それぞれ接続されている。
整数グリッド上の相違度算出モジュール４１、２次曲面係数算出モジュール４２、位置算出モジュール４３、行列Ｂ算出モジュール４４、（５３）式算出モジュール４６は、それぞれ、図２に示した第２−１の実施の形態の整数グリッド上の相違度算出モジュール２１、２次曲面係数算出モジュール２２、位置算出モジュール２３、行列Ｂ算出モジュール２４、（５３）式算出モジュール２５に対応する。
【００７１】
図４に示すように、行列Ｂ算出モジュール４４が行列Ｂを算出した後に、（７１）式正則性検証モジュール４５が前記（７１）式の正則性を検証する。又は、（７１）式正則性検証モジュール４５は、行列Ｂのランクが６以上であることを検証するとしてもよい。
前記（７１）式が正則な場合は、図２に示した第２−１の実施の形態と同様である。
前記（７１）式が正則ではない場合は、（７１）式正則性検証モジュール４５は整数グリッドを増加させる。整数グリッドは前記（７１）式が正則になるまで増加させることになる。整数グリッドを増加させる場合、もし、以前の整数グリッドが１直線上に並んでいる場合には、その直線以外の点を増加させることで、Ｂのランクを増加させることができる。
【００７２】
又は、図５に示すように、第３の実施の形態を変形したものであってもよい。
図５に示す第３の実施の形態は、整数グリッド上の相違度算出モジュール５１、２次曲面係数算出モジュール５２、位置算出モジュール５３、行列Ｂ、Ｗ算出モジュール５４、（７２）式正則性検証モジュール５５、（５４）式算出モジュール５６を有している。整数グリッド上の相違度算出モジュール５１は２次曲面係数算出モジュール５２、行列Ｂ、Ｗ算出モジュール５４、（７２）式正則性検証モジュール５５に、２次曲面係数算出モジュール５２は整数グリッド上の相違度算出モジュール５１、位置算出モジュール５３、（５４）式算出モジュール５６に、位置算出モジュール５３は２次曲面係数算出モジュール５２に、行列Ｂ、Ｗ算出モジュール５４は整数グリッド上の相違度算出モジュール５１、（７２）式正則性検証モジュール５５に、（７２）式正則性検証モジュール５５は整数グリッド上の相違度算出モジュール５１、行列Ｂ、Ｗ算出モジュール５４、（５４）式算出モジュール５６に、（５４）式算出モジュール５６は（７２）式正則性検証モジュール５５、２次曲面係数算出モジュール５２に、それぞれ接続されている。
整数グリッド上の相違度算出モジュール５１、２次曲面係数算出モジュール５２、位置算出モジュール５３、（５４）式算出モジュール５６は、それぞれ、図３に示した第２−２の実施の形態の整数グリッド上の相違度算出モジュール３１、２次曲面係数算出モジュール３２、位置算出モジュール３３、（５４）式算出モジュール３６に対応する。また、行列Ｂ、Ｗ算出モジュール５４は、行列Ｂ算出モジュール３４、重み行列Ｗ算出モジュール３５に対応する。
整数グリッド位置から、行列Ｂ、Ｗ算出モジュール５４は重み行列Ｗも算出し、さらに、（７２）式正則性検証モジュール５５は前記（７２）式の正則性を検証する。
前記（７２）式が正則な場合は、図３に示した第２−２の実施の形態と同様である。
前記（７２）式が正則ではない場合は、整数グリッドを増加させるか、２次曲面近似領域πを広げる。
【００７３】
＜４．７．第４の実施の形態＞
前述の第２の実施の形態の（５３）式算出モジュール２５、第３の実施の形態の（５３）式算出モジュール４６は、毎回（５３）式を計算していたが、前記（５３）式は、整数グリッド位置が決まれば、一意に決定することができるものである。
そこで、代表的な整数グリッド位置に対して、前記（５３）式を予め算出しておくこともできる。
第４の実施の形態は、図６に示すように、整数グリッド上の相違度算出モジュール６１、２次曲面係数算出モジュール６２、位置算出モジュール６３、整数グリッド位置と（５３）式のテーブル６４、行列Ｂ算出モジュール６５、（７１）式正則性検証モジュール６６、（５３）式算出モジュール６７を有している。整数グリッド上の相違度算出モジュール６１は２次曲面係数算出モジュール６２、整数グリッド位置と（５３）式のテーブル６４、（７１）式正則性検証モジュール６６に、２次曲面係数算出モジュール６２は整数グリッド上の相違度算出モジュール６１、整数グリッド位置と（５３）式のテーブル６４、（５３）式算出モジュール６７、位置算出モジュール６３に、位置算出モジュール６３は２次曲面係数算出モジュール６２に、整数グリッド位置と（５３）式のテーブル６４は整数グリッド上の相違度算出モジュール６１、行列Ｂ算出モジュール６５、２次曲面係数算出モジュール６２に、行列Ｂ算出モジュール６５は整数グリッド位置と（５３）式のテーブル６４、（７１）式正則性検証モジュール６６に、（７１）式正則性検証モジュール６６は整数グリッド上の相違度算出モジュール６１、行列Ｂ算出モジュール６５、（５３）式算出モジュール６７に、（５３）式算出モジュール６７は（７１）式正則性検証モジュール６６、２次曲面係数算出モジュール６２に、それぞれ接続されている。
整数グリッド上の相違度算出モジュール６１、２次曲面係数算出モジュール６２、位置算出モジュール６３、行列Ｂ算出モジュール６５、（７１）式正則性検証モジュール６６、（５３）式算出モジュール６７は、それぞれ、図４に示した第３の実施の形態の整数グリッド上の相違度算出モジュール４１、２次曲面係数算出モジュール４２、位置算出モジュール４３、行列Ｂ算出モジュール４４、（７１）式正則性検証モジュール４５、（５３）式算出モジュール４６に対応する。
【００７４】
図６に示すように、整数グリッド位置と（５３）式のテーブル６４に、整数グリッド位置と、前記（５３）式のテーブルを予め用意しておく。
入力整数グリッド位置がテーブル内にある場合には、そのテーブルから前記（５３）式を読み出して、２次曲面を算出すればよい。つまり、２次曲面係数算出モジュール６２が整数グリッド位置と（５３）式のテーブル６４から前記（５３）式の算出結果を受け取る。行列Ｂ算出モジュール６５、（７１）式正則性検証モジュール６６、（５３）式算出モジュール６７による処理は省略される。
入力整数グリッド位置がテーブル内にない場合、図４に示した第３の実施の形態と同様の動作となる。つまり、行列Ｂ算出モジュール６５、（７１）式正則性検証モジュール６６、（５３）式算出モジュール６７による処理が、図４に示した第３の実施の形態と同様に行われる。
【００７５】
又は、図７に示すように、第４の実施の形態を変形したものであってもよい。
図７に示す第４の実施の形態は、整数グリッド上の相違度算出モジュール７１、２次曲面係数算出モジュール７２、位置算出モジュール７３、整数グリッド位置と（５４）式のテーブル７４、行列Ｂ、Ｗ算出モジュール７５、（７２）式正則性検証モジュール７６、（５４）式算出モジュール７７を有している。整数グリッド上の相違度算出モジュール７１は２次曲面係数算出モジュール７２、整数グリッド位置と（５４）式のテーブル７４、（７２）式正則性検証モジュール７６に、２次曲面係数算出モジュール７２は整数グリッド上の相違度算出モジュール７１、位置算出モジュール７３、（５４）式のテーブル７４、（５４）式算出モジュール７７に、位置算出モジュール７３は２次曲面係数算出モジュール７２に、整数グリッド位置と（５４）式のテーブル７４は整数グリッド上の相違度算出モジュール７１、２次曲面係数算出モジュール７２、行列Ｂ、Ｗ算出モジュール７５に、行列Ｂ、Ｗ算出モジュール７５は整数グリッド位置と（５４）式のテーブル７４、（７２）式正則性検証モジュール７６に、（７２）式正則性検証モジュール７６は行列Ｂ、Ｗ算出モジュール７５、整数グリッド上の相違度算出モジュール７１、（５４）式算出モジュール７７に、（５４）式算出モジュール７７は２次曲面係数算出モジュール７２、（７２）式正則性検証モジュール７６にそれぞれ接続されている。
整数グリッド上の相違度算出モジュール７１、２次曲面係数算出モジュール７２、位置算出モジュール７３、行列Ｂ、Ｗ算出モジュール７５、（７２）式正則性検証モジュール７６、（５４）式算出モジュール７７は、それぞれ、図５に示した第３の実施の形態の整数グリッド上の相違度算出モジュール５１、２次曲面係数算出モジュール５２、位置算出モジュール５３、行列Ｂ、Ｗ算出モジュール５４、（７２）式正則性検証モジュール５５、（５４）式算出モジュール５６に対応する。
整数グリッド位置と（５４）式のテーブル７４に、整数グリッド位置と、重み係数行列である前記（５４）式のテーブルを予め用意しておく。整数グリッド位置と（５４）式のテーブル７４は、図６に示した整数グリッド位置と（５３）式のテーブル６４と同様の処理を行う。
つまり、入力整数グリッド位置がテーブル内にある場合には、そのテーブルから前記（５４）式を読み出して、２次曲面を算出すればよい。つまり、２次曲面係数算出モジュール７２が整数グリッド位置と（５４）式のテーブル７４から前記（５４）式の算出結果を受け取る。行列Ｂ、Ｗ算出モジュール７５、（７２）式正則性検証モジュール７６、（５４）式算出モジュール７７による処理は省略される。
入力整数グリッド位置がテーブル内にない場合、図５に示した第３の実施の形態と同様の動作となる。つまり、行列Ｂ、Ｗ算出モジュール７５、（７２）式正則性検証モジュール７６、（５４）式算出モジュール７７による処理が、図５に示した第３の実施の形態と同様に行われる。
【００７６】
さらに、図８に示すように、第４の実施の形態を変形したものであってもよい。
図８に示す第４の実施の形態は、整数グリッド上の相違度算出モジュール８１、２次曲面係数算出モジュール８２、位置算出モジュール８３、整数グリッド位置と（５３）式のテーブル８４を有している。整数グリッド上の相違度算出モジュール８１は２次曲面係数算出モジュール８２、整数グリッド位置と（５３）式のテーブル８４に、２次曲面係数算出モジュール８２は整数グリッド上の相違度算出モジュール８１、整数グリッド位置と（５３）式のテーブル８４、位置算出モジュール８３に、位置算出モジュール８３は２次曲面係数算出モジュール８２に、整数グリッド位置と（５３）式のテーブル８４は整数グリッド上の相違度算出モジュール８１、２次曲面係数算出モジュール８２に、それぞれ接続されている。
整数グリッド上の相違度算出モジュール８１、２次曲面係数算出モジュール８２、位置算出モジュール８３は、それぞれ、図４に示した第３の実施の形態の整数グリッド上の相違度算出モジュール４１、２次曲面係数算出モジュール４２、位置算出モジュール４３に対応する。
整数グリッド位置と（５３）式のテーブル８４は、図６で示した整数グリッド位置と（５３）式のテーブル６４と同様に、整数グリッド位置と前記（５３）式との対応関係を記憶するテーブルである。
図８に示すように、整数グリッド位置と（５３）式のテーブル８４は、テーブル内の整数グリッド位置に合うように、整数グリッドを増加させ、その値を整数グリッド上の相違度算出モジュール８１にフィードバックさせてもよい。テーブルに記述のない余った整数グリッドは捨てて計算すればよい。
そして、整数グリッド位置がテーブル内に存在すれば、対応する前記（５３）式の算出結果を２次曲面係数算出モジュール８２に渡す。
【００７７】
さらに、図９に示すように、第４の実施の形態を変形したものであってもよい。
図９に示す第４の実施の形態は、整数グリッド上の相違度算出モジュール９１、２次曲面係数算出モジュール９２、位置算出モジュール９３、整数グリッド位置と（５４）式のテーブル９４を有している。整数グリッド上の相違度算出モジュール９１は２次曲面係数算出モジュール９２、整数グリッド位置と（５４）式のテーブル９４に、２次曲面係数算出モジュール９２は整数グリッド上の相違度算出モジュール９１、整数グリッド位置と（５４）式のテーブル９４、位置算出モジュール９３に、位置算出モジュール９３は２次曲面係数算出モジュール９２に、整数グリッド位置と（５４）式のテーブル９４は整数グリッド上の相違度算出モジュール９１、２次曲面係数算出モジュール９２に、それぞれ接続されている。
整数グリッド上の相違度算出モジュール９１、２次曲面係数算出モジュール９２、位置算出モジュール９３は、それぞれ、図５に示した第３の実施の形態の整数グリッド上の相違度算出モジュール５１、２次曲面係数算出モジュール５２、位置算出モジュール５３に対応する。
整数グリッド位置と（５４）式のテーブル９４は、図７で示した整数グリッド位置と（５４）式のテーブル７４と同様に、整数グリッド位置と前記（５４）式との対応関係を記憶するテーブルである。整数グリッド位置と（５４）式のテーブル９４は、図７に示した整数グリッド位置と（５４）式のテーブル７４と同様の処理を行う。
つまり、整数グリッド位置と（５４）式のテーブル９４は、テーブル内の整数グリッド位置に合うように、整数グリッドを増加させ、その値を整数グリッド上の相違度算出モジュール９１にフィードバックさせてもよい。テーブルに記述のない余った整数グリッドは捨てて計算すればよい。
そして、整数グリッド位置がテーブル内に存在すれば、対応する前記（５４）式の算出結果を２次曲面係数算出モジュール９２に渡す。
【００７８】
＜４．８．第５の実施の形態＞
＜４．１．３＞で示したように、多次元に拡張することも可能である。これまでの２次曲面係数は、ａ，ｂ，ｃ，ｄ，ｅであったが、多次元拡張する場合には、（７３）式を求めればよい。
【数７３】
2008209971.tif 000088
ただし係数ｆは算出不要である。
さらに、＜４．１．３＞では、行列Ａが半正定値（Ａの固有値が全て０以上）のとき、（７４）式を用いてｈを求めることをしたが、もっと単純にしてもよい。
【数７４】
2008209971.tif 000089
すなわち、行列Ａが正定値（Ａの固有値が全て０より大）のとき、（７４）式を用いてｈを求めるようにしてもよい。
実際の計算上は、逆行列を求めて乗算するよりも、他の高速手法で連立方程式を求めたほうが早いため、（７４）の式を忠実に処理する必要はない。
【００７９】
＜４．９．第６の実施の形態＞
さらに高精度化方式を実装してもよい。
第６の実施の形態は、第１の実施の形態に高精度化方式を実装したものである。
第６の実施の形態は、図１０に示すように、上下、左右方向画素ずらしモジュール１０１、第１の実施の形態の構成１０２、第１の実施の形態の構成１０３、（５１）式算出モジュール１０４を有している。上下、左右方向画素ずらしモジュール１０１は第１の実施の形態の構成１０３に、第１の実施の形態の構成１０２は（５１）式算出モジュール１０４に、第１の実施の形態の構成１０３は上下、左右方向画素ずらしモジュール１０１、（５１）式算出モジュール１０４に、（５１）式算出モジュール１０４は第１の実施の形態の構成１０２、第１の実施の形態の構成１０３に、それぞれ接続されている。
上下、左右方向画素ずらしモジュール１０１は、入力画像２を上下、左右方向にそれぞれ０．５画素ずらした画像を作成する。
この画像と入力画像１を用いて、第１の実施の形態の構成１０３は、第１の実施の形態による処理を行い、ｔｓを求める。
また、第１の実施の形態の構成１０２は、入力画像１と入力画像２とから第１の実施の形態と同様にして、ｔｏを求める。
そして、第１の実施の形態の構成１０２からのｔｏ、第１の実施の形態の構成１０３からのｔｓを用いて、（５１）式算出モジュール１０４が最終的な移動位置を（５１）式で求める。
【数５１】
2008209971.tif 000090
【００８０】
ここでは、変換が２次元の例を示したが、Ｎ次元の場合でも同様である。前記（４５）式を用いて、最終的なｈを求めればよい。
あるいは、前記（４７）式、前記（４８）式、又は前記（４９）式を用いてもよい。これらの式を用いる場合、２次元であれば、上下方向に０．５画素ずらした画像を１種、左右方向に０．５画素ずらした画像を１種作成することになる。図１１では、前記（４６）式を用いる第６の実施の形態の変形例を示している。
第６の実施の形態の変形例は、図１１に示すように、上下方向画素ずらしモジュール１１１、左右方向画素ずらしモジュール１１２、第１の実施の形態の構成１１３、第１の実施の形態の構成１１４、第１の実施の形態の構成１１５、（４６）式算出モジュール１１６を有している。上下方向画素ずらしモジュール１１１は第１の実施の形態の構成１１４に、左右方向画素ずらしモジュール１１２は第１の実施の形態の構成１１５に、第１の実施の形態の構成１１３は（４６）式算出モジュール１１６に、第１の実施の形態の構成１１４は上下方向画素ずらしモジュール１１１、（４６）式算出モジュール１１６に、第１の実施の形態の構成１１５は左右方向画素ずらしモジュール１１２、（４６）式算出モジュール１１６に、（４６）式算出モジュール１１６は第１の実施の形態の構成１１３、第１の実施の形態の構成１１４、第１の実施の形態の構成１１５に、それぞれ接続されている。
上下方向画素ずらしモジュール１１１は、入力画像２を上下方向にそれぞれ０．５画素ずらした画像を作成する。
この画像と入力画像１を用いて、第１の実施の形態の構成１１４は、第１の実施の形態による処理を行い、ｔｓ１を求める。
左右方向画素ずらしモジュール１１２は、入力画像２を左右方向にそれぞれ０．５画素ずらした画像を作成する。
この画像と入力画像１を用いて、第１の実施の形態の構成１１５は、第１の実施の形態による処理を行い、ｔｓ２を求める。
また、第１の実施の形態の構成１１３は、入力画像１と入力画像２から第１の実施の形態と同様にして、ｔｏを求める。
そして、第１の実施の形態の構成１１３からのｔｏ、第１の実施の形態の構成１１４からのｔｓ１、第１の実施の形態の構成１１５からのｔｓ２を用いて、（４６）式算出モジュール１１６が最終的な移動位置を前記（４６）式で求める。
【００８１】
＜４．１０．第７の実施の形態＞
さらに、具体的な応用例を示す。
特許文献３には、画像を重ね合わせて、パノラマ合成を行う技術の例が記載されている。
特許文献３記載の技術では、１３×１３の領域毎に、Ｌｕｃａｓ−Ｋａｎａｄｅ法と呼ばれる画像位置合わせ技術を用いて、例えば図１２に示すような２枚の画像（画像Ａ１２１、画像Ｂ１２２）の位置合わせをサブピクセル単位で行う。このような領域が複数存在していれば、２枚の画像で各領域の中心点の対応を取ることができる。このように対応が取れている点ペアが４点以上あれば、特許文献３記載の数式１６に示されるような射影変換を求めることができる。この射影変換を用いて、２枚の画像（画像Ａ１２１、画像Ｂ１２２）のパノラマ合成（合成画像１２３）をすることができる。
Ｌｕｃａｓ−Ｋａｎａｄｅ法の代わりに、前述の実施の形態の手法を用いることで、同等の処理を行うことができる。
なお、Ｌｕｃａｓ−Ｋａｎａｄｅ法は勾配法の一種である。勾配法は２画像間の２乗誤差を最小化させる特性をもつ手法である。本実施の形態は、２画像間の距離（相違度）に様々な手法を用いることができるため、Ｌｕｃａｓ−Ｋａｎａｄｅ法のような勾配法よりも柔軟で優れていると言える。
【００８２】
＜４．１１．第８の実施の形態＞
別の応用例を示す。動画圧縮への応用例である。
動画圧縮では、局所的な領域（マクロブロック）毎に、隣接するフレーム内で似ている領域を探索し、似ている領域への動きベクトルを符号化することによって、圧縮率を上げる処理が行われる。これは、主に動き補償方式で採用されている。
Ｈ．２６４（ＩＴＵ−Ｔ勧告Ｈ．２６４）等では、上記動きベクトルの単位がサブピクセルとなっている。すなわち、１画素よりも小さな単位で動きベクトルを推定する。本実施の形態を動きベクトル推定に用いることができる。
まず、図１３に示すように、画像Ａ１３１のマクロブロック内の画像と似ている（相違度が小さい）画像Ｂ１３２内の領域を探索する。
最初、この探索は、整数画素単位で行う。
整数画素単位で相違度が最小となる動きベクトルが求まったとする。
次に、その周囲８画素の相違度を算出する。
第１の実施の形態の手法を用いて、サブピクセル単位の動きベクトルを算出することが可能となる。
サブピクセル単位は、１／８等の単位で量子化する必要がある場合には、ｔの各要素をその単位（例えば１／８）で量子化（丸め処理）すればよい。
【００８３】
＜４．１２．実施の形態の補足＞
以上の実施の形態で、中心の整数グリッドが最小の相違度となるように予め移動させておくことを前提としてきたが、これは必ずしも必要ない。
中心の整数グリッドが整数の範囲で最小であれば、精度が高まるというだけであるため、この前提を外してもよい。
行列Ａの行列式が０の場合には、改めて、最小の相違度となる整数グリッドを探してもよい。
【００８４】
＜５．本実施の形態による効果＞
＜５．１．演算量＞
＜５．１．１．グリッド点の量＞
グリッド点の量は、特許文献１に記載の技術で２５点であるのに対し、本実施の形態では９点でも構わない。
ここで、グリッド点の量が全体の処理負荷を決定する。なぜなら、１グリッドの相違度を算出するためには領域Ｗの画素数分程度の積和演算量が必要であるからである。Ｗの画素数は少ない場合でも１００画素程度ある。通常、１０００画素程度である。それに対して移動位置を算出するアルゴリズムに必要な演算負荷は数十回の積和演算である。
よって、グリッド点の量に起因する演算量は、９／２５程度に改善される。しかもこの演算量が支配的なので、全体の演算量も同程度に改善される。
【００８５】
＜５．１．２．最小値算出負荷＞
本実施の形態では、直線毎に最小値を算出する負荷がない。また、ｘ＝±１の直線上でｙ＝±３の範囲に最小値がない場合や、ｙ＝±１の直線上でｘ＝±３の範囲に最小値がない場合でも推定可能である。
【００８６】
＜５．２．グリッド点設定の柔軟性＞
本実施の形態では、グリッド点は、行列Ｂのランクが６以上であれば何でもよい。非常に柔軟に設定可能である。つまり、たまたま存在しているグリッド点の相違度を利用して位置合わせが可能である。例えば、階層的に位置合わせを行う場合、グリッド点の位置が非均等となる可能性が高いが、そのような場合でも利用可能である。
階層的に位置合わせを行う方式の例としては、信学技報（ＩＥ８１−５４）古閑他「会議テレビ信号の動き補償フレーム間符号化」を挙げることができる。
【００８７】
＜５．３．推定精度＞
推定精度に関して、実験を行った結果を示す。
実験画像では、サブピクセル単位で画素位置を人工的にずらして評価を行う。５１２×５１２画素の画像を用いて、これに標準偏差σ＝７のガウシアンフィルタ（入力光学系のＰＳＦを想定している）をかけて帯域制限した後に１０対１にサブサンプリングして１／１０の縮小を行う。こうして作成した縮小画像の中心の３０×３０画素を用いて実験を行った。なお、サブサンプリング位置を１画素ずらすことによって、縮小画像では０．１画素の画素移動を発生させることができる。
実験時には、画素位置移動を−０．５〜０．５の範囲で０．１刻みで発生させる。各手法で位置移動を推定して、推定値のＲＭＳＥを算出する。
実験結果を図１４に示す。横軸に画像の種類を示し、縦軸にＲＭＳＥを示している。各画像に、特許文献１記載の技術、本実施の形態、特許文献１記載の技術を高精度化したもの、本実施の形態の高精度化１、本実施の形態の高精度化２の結果を示している。複数の画像でＲＭＳＥを計測した結果が図１４である。
さらに、各画像のＲＭＳＥをさらに平均したものが図１５に示した表である。
特許文献１記載の技術よりも本実施の形態の精度が高い（ＲＭＳＥが小さい）ことが分かる。本実施の形態の演算負荷は特許文献１記載の技術の演算負荷の９／２５であることを考慮すると、本実施の形態の効果は高い。
さらに、高精度化を行った場合、特許文献１記載の技術と本実施の形態の結果は同じ程度である。演算負荷上は、特許文献１記載の技術では２回補間画像を作成するのに対し、本実施の形態では１回の補間画像の作成で済む。本実施の形態の効果は高い。
さらに高精度化２では、本実施の形態で２回の補間画像を作成した場合の例である。この場合、特許文献１記載の技術よりもさらに高い推定精度を得ることができる。
以上示したように、本実施の形態は特許文献１記載の技術よりも、低演算負荷で、かつ、サブピクセル位置合わせ精度が高いという顕著なメリットがある。
なお、数式を用いて説明したが、数式には、その数式と同等のものを含む。同等のものとは、その数式そのものの他に、最終的な結果に影響を及ぼさない程度の数式の変形、または数式をアルゴリズミックな解法で解くこと等が含まれる。
【００８８】
図２６を参照して、前述の実施の形態のハードウェア構成例について説明する。図２６に示す構成は、例えばパーソナルコンピュータ（ＰＣ）などによって構成されるものであり、スキャナ等のデータ読み取り部２６１７と、プリンタなどのデータ出力部２６１８を備えたハードウェア構成例を示している。
【００８９】
ＣＰＵ（Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）２６０１は、上述の実施の形態において説明した各種のモジュール、すなわち、整数グリッド上の相違度算出モジュール１１、２次曲面係数算出モジュール１２、位置算出モジュール１３等の各モジュールの実行シーケンスを記述したコンピュータ・プログラムにしたがった処理を実行する制御部である。
【００９０】
ＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）２６０２は、ＣＰＵ２６０１が使用するプログラムや演算パラメータ等を格納する。ＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）２６０３は、ＣＰＵ２６０１の実行において使用するプログラムや、その実行において適宜変化するパラメータ等を格納する。これらはＣＰＵバスなどから構成されるホストバス２６０４により相互に接続されている。
【００９１】
ホストバス２６０４は、ブリッジ２６０５を介して、ＰＣＩ（Ｐｅｒｉｐｈｅｒａｌ  Ｃｏｍｐｏｎｅｎｔ  Ｉｎｔｅｒｃｏｎｎｅｃｔ／Ｉｎｔｅｒｆａｃｅ）バスなどの外部バス２６０６に接続されている。
【００９２】
キーボード２６０８、マウス等のポインティングデバイス２６０９は、操作者により操作される入力デバイスである。ディスプレイ２６１０は、液晶表示装置又はＣＲＴ（Ｃａｔｈｏｄｅ  Ｒａｙ  Ｔｕｂｅ）などからなり、各種情報をテキストやイメージ情報として表示する。
【００９３】
ＨＤＤ（Ｈａｒｄ  Ｄｉｓｋ  Ｄｒｉｖｅ）２６１１は、ハードディスクを内蔵し、ハードディスクを駆動し、ＣＰＵ２６０１によって実行するプログラムや情報を記録又は再生させる。ハードディスクは、入力された画像や合成画像などが格納される。さらに、その他の各種のデータ処理プログラム等、各種コンピュータ・プログラムが格納される。
【００９４】
ドライブ２６１２は、装着されている磁気ディスク、光ディスク、光磁気ディスク、又は半導体メモリ等のリムーバブル記録媒体２６１３に記録されているデータ又はプログラムを読み出して、そのデータ又はプログラムを、インタフェース２６０７、外部バス２６０６、ブリッジ２６０５、及びホストバス２６０４を介して接続されているＲＡＭ２６０３に供給する。リムーバブル記録媒体２６１３も、ハードディスクと同様のデータ記録領域として利用可能である。
【００９５】
接続ポート２６１４は、外部接続機器２６１５を接続するポートであり、ＵＳＢ、ＩＥＥＥ１３９４等の接続部を持つ。接続ポート２６１４は、インタフェース２６０７、及び外部バス２６０６、ブリッジ２６０５、ホストバス２６０４等を介してＣＰＵ２６０１等に接続されている。通信部２６１６は、ネットワークに接続され、外部とのデータ通信処理を実行する。データ読み取り部２６１７は、例えばスキャナであり、ドキュメントの読み取り処理を実行する。データ出力部２６１８は、例えばプリンタであり、ドキュメントデータの出力処理を実行する。
【００９６】
なお、図２６に示すハードウェア構成は、１つの構成例を示すものであり、本実施の形態は、図２６に示す構成に限らず、本実施の形態において説明したモジュールを実行可能な構成であればよい。例えば、一部のモジュールを専用のハードウェア（例えば特定用途向け集積回路（Ａｐｐｌｉｃａｔｉｏｎ  Ｓｐｅｃｉｆｉｃ  Ｉｎｔｅｇｒａｔｅｄ  Ｃｉｒｃｕｉｔ：ＡＳＩＣ）等）で構成してもよく、一部のモジュールは外部のシステム内にあり通信回線で接続しているような形態でもよく、さらに図２６に示すシステムが複数互いに通信回線によって接続されていて互いに協調動作するようにしてもよい。また、複写機、ファックス、スキャナ、プリンタ、複合機（多機能複写機とも呼ばれ、スキャナ、プリンタ、複写機、ファックス等の機能を有している）などに組み込まれていてもよい。
【００９７】
なお、説明したプログラムについては、記録媒体に格納して提供してもよく、また、そのプログラムを通信手段によって提供してもよい。その場合、例えば、上記説明したプログラムについて、「プログラムを記録したコンピュータ読み取り可能な記録媒体」の発明として捉えてもよい。
「プログラムを記録したコンピュータ読み取り可能な記録媒体」とは、プログラムのインストール、実行、プログラムの流通などのために用いられる、プログラムが記録されたコンピュータで読み取り可能な記録媒体をいう。
なお、記録媒体としては、例えば、デジタル・バーサタイル・ディスク（ＤＶＤ）であって、ＤＶＤフォーラムで策定された規格である「ＤＶＤ−Ｒ、ＤＶＤ−ＲＷ、ＤＶＤ−ＲＡＭ等」、ＤＶＤ＋ＲＷで策定された規格である「ＤＶＤ＋Ｒ、ＤＶＤ＋ＲＷ等」、コンパクトディスク（ＣＤ）であって、読出し専用メモリ（ＣＤ−ＲＯＭ）、ＣＤレコーダブル（ＣＤ−Ｒ）、ＣＤリライタブル（ＣＤ−ＲＷ）等、光磁気ディスク（ＭＯ）、フレキシブルディスク（ＦＤ）、磁気テープ、ハードディスク、読出し専用メモリ（ＲＯＭ）、電気的消去及び書換可能な読出し専用メモリ（ＥＥＰＲＯＭ）、フラッシュ・メモリ、ランダム・アクセス・メモリ（ＲＡＭ）等が含まれる。
そして、上記のプログラム又はその一部は、上記記録媒体に記録して保存や流通等させてもよい。また、通信によって、例えば、ローカル・エリア・ネットワーク（ＬＡＮ）、メトロポリタン・エリア・ネットワーク（ＭＡＮ）、ワイド・エリア・ネットワーク（ＷＡＮ）、インターネット、イントラネット、エクストラネット等に用いられる有線ネットワーク、あるいは無線通信ネットワーク、さらにこれらの組合せ等の伝送媒体を用いて伝送させてもよく、また、搬送波に乗せて搬送させてもよい。
さらに、上記のプログラムは、他のプログラムの一部分であってもよく、あるいは別個のプログラムと共に記録媒体に記録されていてもよい。また、複数の記録媒体に分割して
記録されていてもよい。また、圧縮や暗号化など、復元可能であればどのような態様で記録されていてもよい。
【図面の簡単な説明】
【００９８】
【図１】第１の実施の形態の構成例についての概念的なモジュール構成図である。
【図２】第２−１の実施の形態の構成例についての概念的なモジュール構成図である。
【図３】第２−２の実施の形態の別の構成例についての概念的なモジュール構成図である。
【図４】第３の実施の形態の構成例についての概念的なモジュール構成図である。
【図５】第３の実施の形態の別の構成例についての概念的なモジュール構成図である。
【図６】第４の実施の形態の構成例についての概念的なモジュール構成図である。
【図７】第４の実施の形態の別の構成例についての概念的なモジュール構成図である。
【図８】第４の実施の形態の別の構成例についての概念的なモジュール構成図である。
【図９】第４の実施の形態の別の構成例についての概念的なモジュール構成図である。
【図１０】第６の実施の形態の構成例についての概念的なモジュール構成図である。
【図１１】第６の実施の形態の別の構成例についての概念的なモジュール構成図である。
【図１２】第７の実施の形態による合成処理例を示す説明図である。
【図１３】第８の実施の形態による処理例を示す説明図である。
【図１４】実験結果の例を示すグラフである。
【図１５】実験結果の例を示す表である。
【図１６】整数グリッドの例を示す説明図である。
【図１７】サブピクセル精度での移動位置の例を示す説明図である。
【図１８】重み係数を設定する例を示す説明図である。
【図１９】ボロノイ分割結果の例を示す説明図である。
【図２０】重み係数を設定する例を示す説明図である。
【図２１】重み係数を設定する例を示す説明図である。
【図２２】整数グリッドの例を示す説明図である。
【図２３】ＳＳＤ値の変化の例を示す説明図である。
【図２４】近似２次曲面の例を示す説明図である。
【図２５】等ＳＳＤ線の例を示す説明図である。
【図２６】本実施の形態を実現するコンピュータのハードウェア構成例を示すブロック図である。
【符号の説明】
【００９９】
１１…整数グリッド上の相違度算出モジュール
１２…２次曲面係数算出モジュール
１３…位置算出モジュール
２１…整数グリッド上の相違度算出モジュール
２２…２次曲面係数算出モジュール
２３…位置算出モジュール
２４…行列Ｂ算出モジュール
２５…（５３）式算出モジュール
３１…整数グリッド上の相違度算出モジュール
３２…２次曲面係数算出モジュール
３３…位置算出モジュール
３４…行列Ｂ算出モジュール
３５…重み行列Ｗ算出モジュール
３６…（５４）式算出モジュール
４１…整数グリッド上の相違度算出モジュール
４２…２次曲面係数算出モジュール
４３…位置算出モジュール
４４…行列Ｂ算出モジュール
４５…（７１）式正則性検証モジュール
４６…（５３）式算出モジュール
５１…整数グリッド上の相違度算出モジュール
５２…２次曲面係数算出モジュール
５３…位置算出モジュール
５４…行列Ｂ、Ｗ算出モジュール
５５…（７２）式正則性検証モジュール
５６…（５４）式算出モジュール
６１…整数グリッド上の相違度算出モジュール
６２…２次曲面係数算出モジュール
６３…位置算出モジュール
６４…整数グリッド位置と（５３）式のテーブル
６５…行列Ｂ算出モジュール
６６…（７１）式正則性検証モジュール
６７…（５３）式算出モジュール
７１…整数グリッド上の相違度算出モジュール
７２…２次曲面係数算出モジュール
７３…位置算出モジュール
７４…整数グリッド位置と（５４）式のテーブル
７５…行列Ｂ、Ｗ算出モジュール
７６…（７２）式正則性検証モジュール
７７…（５４）式算出モジュール
８１…整数グリッド上の相違度算出モジュール
８２…２次曲面係数算出モジュール
８３…位置算出モジュール
８４…整数グリッド位置と（５３）式のテーブル
９１…整数グリッド上の相違度算出モジュール
９２…２次曲面係数算出モジュール
９３…位置算出モジュール
９４…整数グリッド位置と（５４）式のテーブル
１０１…上下、左右方向画素ずらしモジュール
１０２…第１の実施の形態の構成
１０３…第１の実施の形態の構成
１０４…（５１）式算出モジュール
１１１…上下方向画素ずらしモジュール
１１２…左右方向画素ずらしモジュール
１１３…第１の実施の形態の構成
１１４…第１の実施の形態の構成
１１５…第１の実施の形態の構成
１１６…（４６）式算出モジュール
【図１】
2008209971.tif 000112
【図２】
2008209971.tif 000113
【図３】
2008209971.tif 000114
【図４】
2008209971.tif 000115
【図５】
2008209971.tif 000116
【図６】
2008209971.tif 000117
【図７】
2008209971.tif 000118
【図８】
2008209971.tif 000119
【図９】
2008209971.tif 000120
【図１０】
2008209971.tif 000121
【図１１】
2008209971.tif 000122
【図１２】
2008209971.tif 000123
【図１３】
2008209971.tif 000124
【図１４】
2008209971.tif 000125
【図１５】
2008209971.tif 000126
【図１６】
2008209971.tif 000127
【図１７】
2008209971.tif 000128
【図１８】
2008209971.tif 000129
【図１９】
2008209971.tif 000130
【図２０】
2008209971.tif 000131
【図２１】
2008209971.tif 000132
【図２２】
2008209971.tif 000133
【図２３】
2008209971.tif 000134
【図２４】
2008209971.tif 000135
【図２５】
2008209971.tif 000136
【図２６】
2008209971.tif 000137