画像照合装置、画像照合方法及びコンピュータプログラム
【課題】画像照合を高速に実施するとともに、画像の類似性を高精度に判定することが可能な画像照合装置、画像照合方法を提供する。
【解決手段】画像照合装置１は、第１画像及び第２画像を複数のブロックに分割する画像分割部２２、２５と、第１画像及び第２画像について、各ブロックに含まれる画素の正規化画素値の総和及び二乗和に基づき、第１画像の特徴ベクトル及び第２画像の特徴ベクトルを算出する特徴ベクトル算出部２３、２６と、ラグランジュの未定乗算決定法を利用して導いた正規化相互相関の上限値を算出する式を用いて、第１画像の特徴ベクトル及び第２画像の特徴ベクトルに基づき、第１画像と第２画像の正規化相互相関の上限値を算出する上限値算出部２７と、上限値が第１の閾値以上であるか否かに基づいて第１画像と第２画像を照合する第１照合部２８とを有する。
【解決手段】画像照合装置１は、第１画像及び第２画像を複数のブロックに分割する画像分割部２２、２５と、第１画像及び第２画像について、各ブロックに含まれる画素の正規化画素値の総和及び二乗和に基づき、第１画像の特徴ベクトル及び第２画像の特徴ベクトルを算出する特徴ベクトル算出部２３、２６と、ラグランジュの未定乗算決定法を利用して導いた正規化相互相関の上限値を算出する式を用いて、第１画像の特徴ベクトル及び第２画像の特徴ベクトルに基づき、第１画像と第２画像の正規化相互相関の上限値を算出する上限値算出部２７と、上限値が第１の閾値以上であるか否かに基づいて第１画像と第２画像を照合する第１照合部２８とを有する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、画像照合装置、画像照合方法及びコンピュータプログラムに関し、特に、二つの画像が類似するか否かを判定する画像照合装置、画像照合方法及びコンピュータプログラムに関する。
【背景技術】
【０００２】
画像照合技術は、テンプレートマッチング、ブロック動き補償、画像圧縮、ステレオビジョン等の、多くの画像を用いたアプリケーションにおける重要な技術の一つである。最近は、画像の近似コピー検出、画像マイニング及びその画像アノテーションへの応用等、大規模な画像／映像データベースにおける画像照合技術の新しいアプリケーションについても研究が開始されている。
【０００３】
一般に、画像照合の方法として、画素値の差の絶対値の合計（ＳＡＤ：Sum of Absolute Difference）、画素値の差の二乗の合計（ＳＳＤ：Sum of Squared Difference）、正規化相互相関（ＮＣＣ：Normalized Cross-Correlation）等が用いられている。特に、ＮＣＣは、画像間の輝度のずれ、コントラストの違い等による影響が少なく、輝度の補正等がなされた画像間の照合にも適している。例えば、水平方向をｘ軸、垂直方向をｙ軸とする画像Ｉ1と画像Ｉ2についてのＮＣＣ値は、式（１）で定義される。
【数１】
ここで、φx,yは画像Ｉ1の座標（ｘ、ｙ）の画素の画素値であり、ψx,yは画像Ｉ2の座標（ｘ、ｙ）の画素の画素値であり、φaveは画像Ｉ1の全画素の画素値の平均値であり、ψaveは画像Ｉ2の全画素の画素値の平均値である。式（１）に示すように、ＮＣＣ値を算出するためには、対応する画素毎に画素値の相関を計算する必要があるため、特に、照合する画像のサイズが大きい場合には、計算量が大きくなり、算出処理に多大な時間が必要となる。
【０００４】
そこで、例えば、非特許文献１には、ＮＣＣに基づくテンプレートマッチング処理を高速化するための技術が開示されている。非特許文献１に開示された技術は、画像全体に対してテンプレートをずらしながらスキャンする場合、畳み込み演算を用いて各位置での照合結果を求める。一方、空間領域間の畳み込み演算は、周波数領域間での単純な掛け算に変換できるので、空間領域のデータを高速フーリエ変換（ＦＦＴ：Fast Fourie Transform）を用いて周波数領域に変換することにより処理を高速化できる。しかし、この技術はある画像に対してテンプレートをスキャンする場合には効果的であるが、複数の画像に対して照合をする場合は各画像に対してＦＦＴ及び逆ＦＦＴを行う必要があり、かえって処理が遅くなるおそれがある。
【０００５】
一方、非特許文献２には、ＮＣＣの演算処理を効率化するための技術として、ＭＳＥＡ（multilevel successive elimination algorithm）が開示されている。非特許文献２に開示された技術では、式（２）に示すCauchy-Schwarzの不等式に基づいてＮＣＣの上限値が算出される。
【数２】
【０００６】
画像Ｉ1と画像Ｉ2についてブロック内画素数がｍである、ｎ個のブロックに分割した場合、式（１）のＮＣＣ値は、ブロックの番号ｉ（１≦ｉ≦ｎ）及びブロック内の画素の番号ｊ（１≦ｊ≦ｍ）を用いて式（３）で表される。
【数３】
ここで、ｘi,j及びｙi,jは各画素の正規化画素値であり、式（４）で定義される。
【数４】
式（２）、（３）から、ＮＣＣの上限値ＵＢMSEAは、式（５）で表される。
【数５】
ここで、ＸＸi及びＹＹiは式（６）で定義される。
【数６】
ＮＣＣ値はこの上限値ＵＢMSEAを越えることがないので、ＮＣＣ値の代わりに、ＮＣＣ値よりも簡易に算出できる上限値ＵＢMSEAを用いて二つの画像の類似性を評価することができる。
【先行技術文献】
【非特許文献】
【０００７】
【非特許文献１】Tsai, D.-M., Lin, C.-T., 2003. Fast normalized cross correlation for defect detection. Pattern Recognition Letters 24 (15), 2625-2631.
【非特許文献２】Wei, S.-D., Lai, S.-H., 2007. Efficient normalized cross correlation based on adaptive multilevel successive elimination. In: Proceedings of the 8th Asian conference on Computer vision - Volume Part I. ACCV’07. Springer-Verlag, Berlin, Heidelberg, pp. 638-646
【発明の概要】
【発明が解決しようとする課題】
【０００８】
非特許文献２に開示された技術を用いることにより、高速に画像照合を実施することができる。しかしながら、非特許文献２に開示されたＮＣＣの上限値は、ＮＣＣ値より大きくなりすぎる場合がある。その場合、照合する二画像のＮＣＣ値が低いにも関わらず、その二画像の類似性が高いと判定するおそれがあった。
【０００９】
そこで、本発明の目的は、画像照合を高速に実施するとともに、画像の類似性を高精度に判定することが可能な画像照合装置、画像照合方法及びそのような画像照合方法をコンピュータに実行させるコンピュータプログラムを提供することにある。
【課題を解決するための手段】
【００１０】
本発明に係る画像照合装置は、第１画像を複数のブロックに分割し、第２画像を第１画像と同数のブロックに分割する画像分割部と、第１画像及び第２画像のそれぞれについて、各ブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づき、第１画像の特徴ベクトル及び第２画像の特徴ベクトルを算出する特徴ベクトル算出部と、ラグランジュの未定乗算決定法を利用して２つの画像の正規化相互相関値を求める式から導いた２つの画像の正規化相互相関の上限値を算出する式を用いて、算出された第１画像の特徴ベクトル及び第２画像の特徴ベクトルに基づき、第１画像と第２画像の正規化相互相関の上限値を算出する上限値算出部と、上限値が第１の閾値以上であるか否かに基づいて第１画像と第２画像を照合する第１照合部と、を有する。
【００１１】
さらに、本発明に係る画像照合装置において、特徴ベクトル算出部は、ｍをブロック内の画素数、ｎをブロックの数、φi,jを第１画像におけるｉ番目のブロック内のｊ番目の画素の画素値、φaveを第１画像の全画素の画素値の平均値として、次のベクトルξxを第１画像の特徴ベクトルとして算出し、
【数７】
ψi,jを第２画像におけるｉ番目のブロック内のｊ番目の画素の画素値、ψaveを第２画像の全ての画素の画素値の平均値として、次のベクトルξyを第２画像の特徴ベクトルとして算出し、
【数８】
上限値算出部は、第１画像の特徴ベクトルと第２画像の特徴ベクトルの内積を上限値として算出することが好ましい。
【００１２】
さらに、本発明に係る画像照合装置において、上限値が第１の閾値以上である第１画像と第２画像の正規化相互相関値を算出する正規化相互相関値算出部と、正規化相互相関値が第１の閾値以下の値である第２の閾値以上であるか否かに基づいて第１画像と第２画像を詳細に照合する第２照合部と、をさらに有することが好ましい。
【００１３】
また、本発明に係る画像照合方法は、第１画像を複数のブロックに分割するステップと、第１画像について、各ブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づき、第１画像の特徴ベクトルを算出するステップと、第２画像を第１画像と同数のブロックに分割するステップと、第２画像について、各ブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づき、第２画像の特徴ベクトルを算出するステップと、ラグランジュの未定乗算決定法を利用して２つの画像の正規化相互相関値を求める式から導いた２つの画像の正規化相互相関の上限値を算出する式を用いて、算出された第１画像の特徴ベクトル及び第２画像の特徴ベクトルに基づき、第１画像と第２画像の正規化相互相関の上限値を算出するステップと、上限値が第１の閾値以上であるか否かに基づいて第１画像と第２画像を照合するステップと、を含む。
【００１４】
さらに、本発明に係る画像照合方法において、第１画像の特徴ベクトルを算出するステップにおいて、ｍをブロック内の画素数、ｎをブロックの数、φi,jを第１画像におけるｉ番目のブロック内のｊ番目の画素の画素値、φaveを第１画像の全画素の画素値の平均値として、次のベクトルξxを第１画像の特徴ベクトルとして算出し、
【数９】
第２画像の特徴ベクトルを算出するステップにおいて、ｍをブロック内の画素数、ｎをブロックの数、ψi,jを第２画像におけるｉ番目のブロック内のｊ番目の画素の画素値、ψaveを第２画像の全ての画素の画素値の平均値として、次のベクトルξyを第２画像の特徴ベクトルとして算出し、
【数１０】
上限値を算出するステップにおいて、第１画像の特徴ベクトルと第２画像の特徴ベクトルの内積を上限値として算出することが好ましい。
【００１５】
さらに、本発明に係る画像照合方法において、上限値が第１の閾値以上である第１画像と第２画像の正規化相互相関値を算出するステップと、正規化相互相関値が第１の閾値以下の値である第２の閾値以上であるか否かに基づいて第１画像と第２画像を詳細に照合するステップと、をさらに有することが好ましい。
【００１６】
また、本発明に係るコンピュータプログラムは、第１画像を複数のブロックに分割するステップと、第１画像について、各ブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づき、第１画像の特徴ベクトルを算出するステップと、第２画像を第１画像と同数のブロックに分割するステップと、第２画像について、各ブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づき、第２画像の特徴ベクトルを算出するステップと、ラグランジュの未定乗算決定法を利用して２つの画像の正規化相互相関値を求める式から導いた２つの画像の正規化相互相関の上限値を算出する式を用いて、算出された第１画像の特徴ベクトル及び第２画像の特徴ベクトルに基づき、第１画像と第２画像の正規化相互相関の上限値を算出するステップと、上限値が第１の閾値以上であるか否かに基づいて第１画像と第２画像を照合するステップと、をコンピュータに実行させる。
【発明の効果】
【００１７】
本発明によれば、画像照合を高速に実施するとともに、画像の類似性を高精度に判定することが可能な画像照合装置、画像照合方法及びそのような画像照合方法をコンピュータに実行させるコンピュータプログラムを提供することができる。
【図面の簡単な説明】
【００１８】
【図１】本発明を適用した画像照合装置の概略構成図である。
【図２】画像照合処理の動作を示すフローチャートである。
【図３】ブロック毎に分割した画像の例である。
【図４】制御部の他の例を示す概略構成図である。
【図５】詳細な画像照合処理の動作を示すフローチャートである。
【図６】制御部の他の例を示す概略構成図である。
【図７】画像の取得処理の動作を示すフローチャートである。
【図８】画像の照合処理の動作を示すフローチャートである。
【図９】特徴量次元数毎のプレシジョンを表すグラフである。
【図１０】ＤＣＴを用いた低次元特徴量を説明するための概略図である。
【図１１】（ａ）、（ｂ）は、リコールとプレシジョンの関係のグラフである。
【図１２】ＮＩＨを用いた低次元特徴量を説明するための概略図である。
【図１３】（ａ）、（ｂ）は、リコールとプレシジョンの関係のグラフである。
【図１４】（ａ）〜（ｄ）は、特徴量次元数毎のプレシジョンを表すグラフである。
【図１５】（ａ）〜（ｄ）は、特徴量次元数毎のプレシジョンを表すグラフである。
【図１６】特徴量次元数毎の照合処理時間を示す表である。
【図１７】画像照合装置の他の例を示す概略構成図である。
【図１８】画像の符号化処理の動作を示すフローチャートである。
【図１９】インテグラル画像を説明するための概略図である。
【発明を実施するための形態】
【００１９】
以下、本発明に係る画像照合装置、画像照合方法及びコンピュータプログラムについて図を参照しつつ説明する。但し、本発明の技術的範囲はそれらの実施の形態に限定されず、特許請求の範囲に記載された発明とその均等物に及ぶ点に留意されたい。
【００２０】
図１は、本発明を適用した画像照合装置の概略構成を示す図である。図１に示すように、画像照合装置１は、インターフェース部１１、記憶部１２、表示部１３及び制御部２０を有する。以下、画像照合装置１の各部について詳細に説明する。
【００２１】
インターフェース部１１は、例えばインターネット、電話回線網（携帯端末回線網、一般電話回線網を含む）、イントラネット等のネットワークを介して他のコンピュータ等と画像データ及び各種のデータを送受信する通信インターフェースであり、接続するネットワークの通信インターフェース回路を有する。また、インターフェース部１１は、例えばＵＳＢ等のシリアルバスに準じるインターフェース回路を有し、フラッシュメモリ等を接続し、そのフラッシュメモリ等から画像データ及び各種のデータを取得するようにしてもよい。インターフェース部１１は、制御部２０と接続されており、制御部２０により制御される。
【００２２】
記憶部１２は、ＲＡＭ、ＲＯＭ等のメモリ装置、ハードディスク等の固定ディスク装置、又はフレキシブルディスク、光ディスク等の可搬用の記憶装置等を有する。また、記憶部１２には、画像照合装置１の各種処理に用いられるコンピュータプログラム、データベース、テーブル等が格納される。記憶部１２は、制御部２０と接続され、インターフェース部１１を介して取得した画像データを格納するとともに、制御部２０により画像データについてなされた演算結果を格納する。
【００２３】
制御部２０は、複数の画像について特徴ベクトルを算出し、各画像が類似するか否かを判定する。そのために、制御部２０は、第１画像取得部２１、第１画像分割部２２、第１特徴ベクトル算出部２３、第２画像取得部２４、第２画像分割部２５、第２特徴ベクトル算出部２６、上限値算出部２７及び第１照合部２８を有する。また、制御部２０は、インターフェース部１１、記憶部１２及び表示部１３と接続され、インターフェース部１１のデータ送受信制御、記憶部１２の制御、表示部１３の表示制御等を行う。制御部２０は、予め記憶部１２に記憶されているプログラムに基づいて動作する。あるいは、制御部２０は、集積回路、マイクロプロセッサ、ファームウェア等で構成されてもよい。
【００２４】
図２は、画像照合装置１による画像照合処理の動作を示すフローチャートである。以下、図２に示したフローチャートを参照しつつ、画像照合処理の動作を説明する。なお、以下に説明する動作のフローは、予め記憶部１２に記憶されているプログラムに基づき主に制御部２０により画像照合装置１の各要素と協同して実行される。
【００２５】
最初に、第１画像取得部２１は、インターフェース部１１を介して、外部のコンピュータ、フラッシュメモリ等から画像を取得し（以下、第１画像取得部２１が取得した画像を第１画像と称する）、記憶部１２に保存する（ステップＳ２０１）。
【００２６】
次に、第１画像分割部２２は、記憶部１２に保存された第１画像を読み出し、第１画像を所定サイズの複数のブロックに分割する（ステップＳ２０２）。例えば、第１画像分割部２２は、３５２×２４０画素の画像を２２×１５画素の２５６ブロックに分割する。
【００２７】
図３は、ブロック毎に分割した画像の例を示す模式図である。図３の例では、画像３００を、各ブロックの画素数がｍであるｎ個のブロックに分割した画像３１０を示す。
【００２８】
次に、第１特徴ベクトル算出部２３は、第１画像について、分割したブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づく特徴ベクトルを算出し、第１画像と関連付けて記憶部１２に保存する（ステップＳ２０３）。
【００２９】
画像ｘ、画像ｙについてブロック内画素数がｍである、ｎ個のブロックに分割した場合のＮＣＣ値を、ブロックの番号ｉ（１≦ｉ≦ｎ）及びブロック内の画素の番号ｊ（１≦ｊ≦ｍ）を用いて式（１１）で定義する。
【数１１】
ここで、φi,j、ψi,jは、各画素の画素値であり、φaveは画像ｘの全画素の画素値の平均値であり、ψaveは画像ｙの全画素の画素値の平均値である。
【００３０】
画像ｘ及び画像ｙの各画素の正規化画素値ｘi,j及びｙi,jを式（１２）で定義する。
【数１２】
なお、正規化画素値は、厳密には式（１３）とすべきであるが、本実施形態では単純化のため正規化画素値を式（１２）で定義する。
【数１３】
【００３１】
式（１１）のＮＣＣ値は、式（１２）の正規化画素値ｘi,j及びｙi,jを用いて式（１４）で表される。
【数１４】
【００３２】
また、画像ｘのｉ番目のグループにおける正規化画素値ｘi,jの総和及び二乗和並びに画像ｙのｉ番目のグループにおける正規化画素値ｙi,jの総和及び二乗和を式（１５）〜（１８）で定義する。
【数１５】
【数１６】
【数１７】
【数１８】
ここで、正規化画素値の性質上、式（１９）、（２０）が成り立つ。
【数１９】
【数２０】
【００３３】
以下、ラグランジュの未定乗数決定法を利用して式（１４）から２つの画像のＮＣＣの厳密な上限値を算出する式を求める。本実施形態では、式（２１）の関数の極値問題を考える。
【数２１】
ここで、λ1,i、λ2,i、λ3,i、λ4,iは、４ｎ個のラグランジュ未定乗数である。式（２１）において、微分係数ｄΛ＝０（全てのｉ、ｊについて∂Λ/∂ｘi,j＝０、∂Λ/∂yi,j＝０）とすることにより、式（２２）、（２３）の２ｍｎ個の式が得られる。
【数２２】
【数２３】
式（２２）、（２３）から式（２４）が得られる。
【数２４】
【００３４】
ここで、式（２４）が成立すると仮定すると、ｘi,j、ｙi,jは式（２５）で求められる。
【数２５】
この場合、式（２５）の右辺はｉ（ブロックの番号）のみにより定まるので、ｘi,j、ｙi,jは全てのｊ（ブロック内の画素）について同一の値をとることになり、条件が強すぎることになる。一方、式（２４）が成立しないと仮定すると、式（２４）の左辺の行列の行列式から式（２６）が得られる。
【数２６】
また、式（２２）、（２３）、（２６）から、式（２７）が得られる。
【数２７】
この場合、式（２２）と式（２３）は、同じ意味を表すことになる。従って、本実施形態では、式（２２）のみを用いてＮＣＣの上限値を算出する式を求める。
【００３５】
まず、式（２２）から式（２８）が得られる。
【数２８】
式（２８）を式（１７）に代入することにより、式（２９）が得られる。
【数２９】
一方、式（２８）を式（１８）に代入することにより、式（３０）が得られる。
【数３０】
式（２９）より、式（３１）が得られる。
【数３１】
式（３１）を式（３０）に代入することにより、式（３２）が得られる。
【数３２】
式（２８）、（３１）、（３２）を式（１４）に代入することにより、式（３３）が得られる。
【数３３】
【００３６】
つまり、ＮＣＣの上限値ＵＢは式（３４）で与えられ、ＮＣＣの下限値ＬＢは式（３５）で与えられる。
【数３４】
【数３５】
ＮＣＣ値はこの上限値ＵＢを越えることがないので、ＮＣＣ値の代わりに上限値ＵＢを用いて画像ｘと画像ｙの類似性を評価することができる。
【００３７】
一方、式（３４）の上限値ＵＢは、式（３６）に示されるように、式（３７）、（３８）で示される特徴ベクトルξxとξyの内積となる。
【数３６】
【数３７】
【数３８】
式（３７）、（３８）に示されるように、特徴ベクトルξx及びξyは、それぞれ画像ｘ、画像ｙのみから算出することができる。つまり、画像照合装置１は、予め画像ｘ、画像ｙをそれぞれ２ｎ個の要素をもつ特徴ベクトルξx及びξyに変換しておくことにより、その内積を算出して上限値ＵＢを求めることができ、特に、ある画像を複数の画像と照合するときに高速に照合することができる。そこで、第１特徴ベクトル算出部２３は、第１画像の特徴ベクトルとして第１画像について式（３７）の特徴ベクトルξxを算出する。
【００３８】
次に、第２画像取得部２４は、インターフェース部１１を介して、外部のコンピュータ、フラッシュメモリ等から画像を取得し（以下、第２画像取得部２４が取得した画像を第２画像と称する）、記憶部１２に保存する（ステップＳ２０４）。
【００３９】
次に、第２画像分割部２５は、記憶部１２に保存された第２画像を読み出し、第１画像を分割したブロックとそれぞれ同サイズ、同数の複数のブロックに分割する（ステップＳ２０５）。
【００４０】
次に、第２特徴ベクトル算出部２６は、第２画像について、分割したブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づく特徴ベクトル、つまり式（３８）の特徴ベクトルξyを算出し、第２画像と関連付けて記憶部１２に保存する（ステップＳ２０６）。
【００４１】
次に、上限値算出部２７は、式（３６）に示すように特徴ベクトルξx、ξyの内積を算出して、上限値ＵＢを求める（ステップＳ２０７）。
【００４２】
次に、第１照合部２８は、上限値ＵＢが所定の閾値θ1以上であるか否かに基づいて第１画像と前記第２画像を照合する。まず、第１照合部２８は、上限値ＵＢが閾値θ1以上であるか否かを判定する（ステップＳ２０８）。この閾値θ1は、ＮＣＣを用いて二つの画像が類似するか否かを判定するときに、ＮＣＣ値と比較する閾値θ0（−１≪θ0＜１）と同じ値である。あるいは、閾値θ1は、閾値θ0より大きい値としてもよい。
【００４３】
第１照合部２８は、上限値ＵＢが閾値θ1以上である場合、第１画像が第２画像と類似すると判定し（ステップＳ２０９）、一連のステップを終了する。一方、第１照合部２８は、上限値ＵＢが閾値θ1未満である場合、第１画像が第２画像と類似しないと判定し（ステップＳ２１０）、一連のステップを終了する。
【００４４】
例えば、画像照合装置１は、類似すると判定した画像のペアを利用者が確認できるように、表示部１３に表示する。あるいは、画像照合装置１は、類似すると判定した画像のペアの情報をインターフェース部１１を介して外部のコンピュータ（不図示）に通知してもよい。
【００４５】
なお、上限値ＵＢと下限値ＬＢの間の範囲が狭いほど、上限値ＵＢは、ＮＣＣ値をより正確に近似していることになる。上限値ＵＢと下限値ＬＢの間の範囲は、式（３４）、（３５）から、式（３９）により求めることができる。
【数３９】
一方、
【数４０】
であり、ここで、Ｇj,k＝ｘi,jｘi,j−ｘi,jｘi,kとすると、式（４１）が成り立つ。
【数４１】
従って、式（４０）は、式（４２）で表される。
【数４２】
ここで、Ｅkはｋにわたる期待値であり、Ｖａｒi(ｘi,j)はｉ番目のブロック内のｘi,jの分散である。同様に、
【数４３】
が成立し、式（３９）は、式（４４）で表される。
【数４４】
【００４６】
つまり、上限値ＵＢと下限値ＬＢの間の範囲を狭くするためには、各ブロック内部の正規化画素値の分散を小さくするか、一ブロック内の画素数ｍを小さくする必要がある。しかし、一ブロック内の画素数ｍを小さくすると、ブロック数ｎが大きくなり、特徴ベクトルξxとξyの次元数が大きくなるため、現実的でない。従って、同一ブロック内の画素の間の正規化画素値の差が可能な限り小さくなるように各ブロックを構成し、各ブロック内部の正規化画素値の分散を小さくすることが好ましい。一般に、画像のマルコフ性により、近傍の画素の（正規化）画素値は、非常に近い値をもつことが知られている。従って、例えばブロックの形状を、長方形、特に正方形に近い形にすることにより、各ブロック内の画素の画素値を近付け、上限値ＵＢをよりＮＣＣ値に近似させることができる。
【００４７】
なお、ＭＳＥＡにより算出されるＮＣＣの上限値ＵＢMSEAは、式（５）から式（４５）のように表される。
【数４５】
本実施形態により算出されるＮＣＣの上限値ＵＢと、ＭＳＥＡにより算出されるＮＣＣの上限値ＵＢMSEAについて、式（３４）と式（４５）のルート内の式を比較することにより式（４６）が成立する。
【数４６】
従って、式（４７）が成立し、本実施形態により算出されるＮＣＣの上限値ＵＢは、ＭＳＥＡにより算出されるＮＣＣの上限値ＵＢMSEA以下となる。
【数４７】
つまり、本実施形態により算出されるＮＣＣの上限値ＵＢは、ＭＳＥＡにより算出されるＮＣＣの上限値ＵＢMSEAより正確にＮＣＣ値を近似していることが理論的に示される。
【００４８】
以上詳述したように、図２に示したフローチャートに従って動作することによって、画像照合装置１は、元の情報量（数万程度の画素数）から大幅に圧縮した（数十〜百程度の）特徴ベクトル（すなわち低次元特徴量）を用いて、画像の類似性を高速かつ高精度に判定することが可能となった。
【００４９】
なお、画像照合装置１は、照合する二画像のサイズが異なる場合、二画像のサイズを同一サイズに正規化することにより、サイズの異なる二画像についても低次元特徴量を用いた画像照合をすることができる。
【００５０】
図４は、制御部の他の例を示す概略構成図である。図４に示す制御部３０は、図１に示す制御部２０の代わりに用いることが可能である。図４に示す制御部３０は、図１に示す制御部２０の各部に加えて、ＮＣＣ値算出部３１と第２照合部３２を有する。
【００５１】
図５は、図４に示す制御部３０を用いる画像照合装置１による詳細な画像照合処理の動作を示すフローチャートである。以下、図５に示したフローチャートを参照しつつ、詳細な画像照合処理の動作を説明する。なお、以下に説明する動作のフローは、予め記憶部１２に記憶されているプログラムに基づき主に制御部３０により画像照合装置１の各要素と協同して実行される。また、図５に示したフローチャートは、図２に示したフローチャートにより、類似すると判定された画像のペアに対して実行される。
【００５２】
最初に、ＮＣＣ値算出部３１は、類似すると判定された第１画像と第２画像について式（１１）を用いてＮＣＣ値を算出する（ステップＳ５０１）。
【００５３】
次に、第２照合部３２は、ＮＣＣ値算出部３１により算出されたＮＣＣ値が所定の閾値θ0（例えば０．９）以上であるか否かに基づいて第１画像と第２画像を詳細に照合する。まず、第２照合部３２は、ＮＣＣ値が閾値θ0以上であるか否かを判定する（ステップＳ５０２）。なお、上述した通り、この閾値θ0は閾値θ1以下の値の閾値である。
【００５４】
第２照合部３２は、ＮＣＣ値が閾値θ0以上である場合、第１画像と第２画像は類似すると判定し（ステップＳ５０３）、一連のステップを終了する。一方、第２照合部２３は、ＮＣＣ値が閾値θ0未満である場合、第１画像と第２画像は類似しないと判定し（ステップＳ５０４）、一連のステップを終了する。
【００５５】
例えば、画像照合装置１は、図５に示したフローチャートにより類似すると判定した画像のペアを表示部１３に表示する。あるいは、画像照合装置１は、図５に示したフローチャートにより類似すると判定した画像のペアの情報をインターフェース部１１を介して外部のコンピュータ（不図示）に通知してもよい。
【００５６】
なお、制御部２０は、第１画像取得部２１と第２画像取得部２４の代わりに、第１画像及び第２画像の何れをも取得できる画像取得部を備えることもできる。また、第１画像分割部２２と第２画像分割部２５の代わりに第１画像及び第２画像の何れをも所定のブロックに分割できる画像分割部を備えることもできる。さらに、第１特徴ベクトル算出部２３と第２特徴ベクトル算出部２６の代わりに第１画像の特徴ベクトル及び第２画像の特徴ベクトルの何れをも算出できる特徴ベクトル算出部を備えることもできる。
【００５７】
また、ステップＳ２０４〜Ｓ２０６の第２画像に対する処理は、ステップＳ２０１〜Ｓ２０３の第１画像に対する処理より前に実施してもよい。また、画像照合装置１が複数のＣＰＵを備えること等により、並列処理が可能な場合には、ステップＳ２０４〜Ｓ２０６の第２画像に対する処理と、ステップＳ２０１〜Ｓ２０３の第１画像に対する処理を並列に実施してもよい。
【００５８】
以上詳述したように、図５に示したフローチャートに従って動作することによって、画像照合装置１は、低次元空間での類似性が高い（ＮＣＣの上限値ＵＢが閾値θ1以上の）画像のペアのうち、オリジナル空間での類似性が低い（ＮＣＣ値が閾値θ0未満の）画像のペアを取り除く（クレンジングする）ことができる。これにより、画像照合装置１は、ある画像のペアに対して低負荷で類似性を判定しておき、その判定において類似性が高いと判定した画像のペアのみについて高負荷で高精度に類似性を判定するので、効率よく画像を照合できるようになった。
【００５９】
図６は、画像照合装置を画像検索サーバに適用する場合の制御部の概略構成図である。図６に示す制御部４０は、図１に示す制御部２０又は図４に示す制御部３０の代わりに用いることが可能である。図６に示す制御部４０は、図４に示す制御部３０の各部に加えて、インデックス判定部４１を有する。なお、以下では、第１画像を事前に記憶部１２に格納された照合元の画像とし、第２画像を外部のコンピュータ（不図示）等から検索を要求された問合せ画像として説明する。
【００６０】
図７は、図６に示す制御部４０を用いる画像照合装置１による画像の取得処理の動作を示すフローチャートである。以下、図７に示したフローチャートを参照しつつ、画像取得処理の動作を説明する。なお、以下に説明する動作のフローは、予め記憶部１２に記憶されているプログラムに基づき主に制御部４０により画像照合装置１の各要素と協同して実行される。
【００６１】
図７に示すフローチャートでは、画像照合装置１は、第１画像の取得のみを実施し、第１画像と第２画像の照合処理は、後述するフローチャートで説明する。図７に示すステップＳ７０１〜Ｓ７０３の処理は、図２に示すステップＳ２０１〜Ｓ２０３の処理と同じであるため、説明を省略する。
【００６２】
第１特徴ベクトル算出部２３は、第１画像について特徴ベクトルを算出すると、算出した特徴ベクトルに基づいて多次元インデックスを作成する（ステップＳ７０４）。
【００６３】
式（４８）、（４９）に示すように、式（３７）、（３８）で求められる特徴ベクトルξx、ξyのノルムは１になる。
【数４８】
【数４９】
つまり、式（５０）に示すように、特徴ベクトルξxとξyの内積は、特徴ベクトルξxとξyのユークリッド距離に変換することができる。
【数５０】
つまり、特徴ベクトルξxとξyのユークリッド距離が式（５１）の条件を満たす場合、特徴ベクトルξxとξyの内積は閾値θ1以上となり、ＮＣＣの上限値ＵＢは閾値θ1以上となる。
【数５１】
【００６４】
従って、画像照合装置１は、特徴ベクトルξxとξyの内積を算出する代わりに、特徴ベクトルξxとξyのユークリッド空間における距離探索を行うことにより、式（５２）の条件を満たす画像のペアを抽出することができる。
【数５２】
【００６５】
この場合、画像照合装置１は、ユークリッド距離を用いた距離探索をサポートする、多次元データに対するインデキシング技術（Bohm, C., Berchtold, S., Keim, D. A., September 2001. Searching in high-dimensional spaces: Index structures for improving the performance of multimedia databases. ACM Computing Surveys 33 (3), 322-373.）を適用することができ、画像照合のさらなる高速化を図ることができる。画像照合装置１は、多次元インデキシング技術として、例えば、木構造を用いるＡＮＮ（Approximate Nearest Neighbors）（Arya, S., Mount, D. M., Netanyahu, N. S., Silverman, R., Wu, A. Y., November 1998. An optimal algorithm for approximate nearest neighbor searching fixed dimensions. J. ACM 45, 891-923.）、ハッシュを用いるＬＳＨ（Locality Sensitive Hashing）（Andoni, A., Indyk, P., 2008. Near-optimal hashing algorithms for approximate nearest neighbor in high dimensions. Commun. ACM 51 (1), 117-122.）、階層的ベクトル量子化（Nister, D., Stewenius, H., 2006. Scalable recognition with a vocabulary tree. In: Proc. of CVPR. Vol. 2. pp. 2161-2168.）、ベクトル量子化及びスカラー量子化（ハミングエンベディング）（Douze, M., Je´gou, H., Sandhawalia, H., Amsaleg, L., Schmid, C., 2009. Evaluation of gist descriptors for web-scale image search. In: Proceeding of the ACM International Conference on Image and Video Retrieval. CIVR ’09. ACM, New York, NY, USA, pp. 19:1-19:8.）等の、高次元空間におけるｋ近傍法（ｋ−ＮＮ探索）を用いた技術を利用できる。
【００６６】
例えば、多次元インデキシング技術としてＬＳＨを利用する場合、第１特徴ベクトル算出部２３は、特徴ベクトルに対して、予め定められた複数のハッシュ関数を用いて複数のハッシュ値を算出し、多次元インデックスとする。
【００６７】
あるいは、多次元インデキシング技術としてＡＮＮを利用する場合、画像照合装置１は、特徴ベクトルについての特徴空間を予め所定の領域に分割しておく。そして、第１特徴ベクトル算出部２３は、特徴ベクトルが属する特徴空間の領域を多次元インデックスとする。
【００６８】
第１特徴ベクトル算出部２３は、多次元インデックスを作成すると、第１画像と関連付けて記憶部１２に格納し（ステップＳ７０５）、一連のステップを終了する。
【００６９】
図８は、画像照合装置１が画像検索サーバとして動作する場合の画像照合処理の動作を示すフローチャートである。以下、図８に示したフローチャートを参照しつつ、画像照合処理の動作を説明する。なお、以下に説明する動作のフローは、予め記憶部１２に記憶されているプログラムに基づき主に制御部４０により画像照合装置１の各要素と協同して実行される。
【００７０】
なお、図８に示すステップＳ８０１〜Ｓ８０３の処理は、図２に示すステップＳ２０４〜Ｓ２０６の処理と同じであるため、説明を省略する。ただし、多次元インデキシング技術としてＬＳＨを利用する場合、第２特徴ベクトル算出部２６は、ステップＳ８０３で第２画像について特徴ベクトルを算出した後、その特徴ベクトルについて、第１画像の特徴ベクトルに対してハッシュ値を算出したのと同一の、複数のハッシュ関数を用いて複数のハッシュ値を算出する。
【００７１】
ステップＳ８０４〜Ｓ８０８の処理は、記憶部１２に格納された全ての第１画像に対して各画像毎に実施される。
【００７２】
ステップＳ８０４において、インデックス判定部４１は、第２特徴ベクトル算出部２６により第２画像について特徴ベクトルが算出されると、記憶部１２から第１画像の特徴ベクトル及び多次元インデックスを読み出す（ステップＳ８０４）。
【００７３】
次に、インデックス判定部４１は、第１画像の特徴ベクトルと第２画像の特徴ベクトルが多次元インデックスに基づく条件を満たしているか否かを判定する（ステップＳ８０５）。
【００７４】
例えば、多次元インデキシング技術としてＬＳＨを利用する場合、インデックス判定部４１は、記憶部１２から読み出した第１画像の特徴ベクトルのハッシュ値が、第２画像の特徴ベクトルの対応するハッシュ値と所定数（例えば、全ハッシュ値のうちの半数）以上一致するか否かにより、多次元インデックスに基づく条件を満たしているか否かを判定する。
【００７５】
あるいは、多次元インデキシング技術としてＡＮＮを利用する場合、インデックス判定部４１は、第１画像の特徴ベクトルが属する特徴空間の領域が、第２画像の特徴ベクトルとのユークリッド距離が（２−２θ1）以下となる部分を含むか否かにより、多次元インデックスに基づく条件を満たしているか否かを判定する。
【００７６】
インデックス判定部４１が多次元インデックスに基づく条件を満たしていると判定した場合、第１照合部２８は、第１画像の特徴ベクトルと第２画像の特徴ベクトルのユークリッド距離が式（５１）の条件を満たすか否かを判定する（ステップＳ８０６）。
【００７７】
第１照合部２８は、式（５１）の条件が満たされている場合、その第１画像を第２画像に類似する画像の候補（以下、類似画像候補と称する）として抽出する（ステップＳ８０７）。
【００７８】
一方、ステップＳ８０５でインデックス判定部４１が多次元インデックスに基づく条件が満たされていないと判定した場合、ステップＳ８０６で第１照合部２８が式（５１）の条件が満たされていないと判定した場合、又は、ステップＳ８０７で第１照合部２８が類似画像候補を抽出した場合、制御部４０は、第２画像を全ての第１画像と比較したか否かを判定する（ステップＳ８０８）。
【００７９】
制御部４０は、まだ第２画像と比較していない第１画像がある場合、ステップＳ８０４〜Ｓ８０７の処理を繰り返し、全ての第２画像と第１画像との比較が完了すると、一連のステップを終了する。
【００８０】
制御部４０は、図８のフローチャートを終了した後、さらに、図５に示したフローチャートを実施し、第２画像の類似画像候補として抽出した第１画像についてＮＣＣ値を用いて詳細に画像照合を行う。そして、制御部４０は、第２画像と類似すると判定した第１画像の情報をインターフェース部１１を介して検索を要求したコンピュータに通知する。あるいは、制御部４０は、図８のフローチャートにより類似画像候補として抽出した第１画像の情報をインターフェース部１１を介して検索を要求したコンピュータに通知してもよい。
【００８１】
以上詳述したように、図８に示したフローチャートに従って動作することによって、画像照合装置１は、画像検索サーバとして動作し、複数の画像に対する画像照合を高速かつ高精度に実施することが可能となった。また、所定の画像を複数の画像と照合する場合に多次元インデックス構造を利用して照合処理を高速化することが可能となった。
【００８２】
以下、本実施形態の低次元変換方法による類似画像の抽出精度について説明する。複数の画像Ｉi1を有する画像セットＩＳ1（ＩＳ1=｛Ｉi1｝）と複数の画像Ｉi2を有する画像セットＩＳ2（ＩＳ2=｛Ｉi2｝）において、ＮＣＣ値が閾値θ0より大きくなる画像のペアの集合ＩＩＰＳは、式（５３）のように表される。
【数５３】
この集合ＩＩＰＳは、本実施形態で算出されるＮＣＣの上限値ＵＢが閾値θ0より大きくなる、つまり式（５４）で表される画像のペアの集合ＩＩＰＳ'に近似できる。
【数５４】
【００８３】
上述した通り、ＮＣＣ値は上限値ＵＢ以下となるので、上限値ＵＢがθ0以下である場合、ＮＣＣ値は確実にθ0以下となり、対応する画像のペアは確実にＩＩＰＳに含まれないと判断できる。つまり、式（５５）の関係が成立し、ＩＩＰＳに含まれる画像のペアがＩＩＰＳ'から取りこぼされることはない。
【数５５】
【００８４】
一方、ＩＩＰＳに対してＩＩＰＳ'が大きくなりすぎると、類似しない画像まで類似画像として抽出していることになり、類似画像の抽出精度が低いということになる。このように、低次元変換の効果及び抽出した画像の類似性は、類似画像の抽出精度（以下、プレシジョンと称する）と、類似画像の非取りこぼし率（以下、リコールと称する）を用いて評価することができる。式（５６）はプレシジョンの算出式であり、式（５７）はリコールの算出式である。
【数５６】
【数５７】
例えば、ＩＩＰＳ'を非常に大きくする（極端な例ではＩＩＰＳ'をＩＩＰＳの全集合とする）ことにより、リコールを容易に１００％にすることができるが、その場合、プレシジョンが低くなる。低次元変換においてはリコールを１００％に保ちつつ、プレシジョンを可能な限り高くすることが好ましい。
【００８５】
図９は、本実施形態及び他の低次元変換方法においてリコールを１００％にしたときのプレシジョンを表すグラフである。図９では、他の低次元変換方法としてＭＳＥＡと、離散コサイン変換（ＤＣＴ：Discrete Cosine Transform）が用いられる。
【００８６】
図９に示すグラフ９００において、グラフ９０１は本実施形態によるプレシジョンを示し、グラフ９０２はＭＳＥＡによるプレシジョンを示し、グラフ９０３はＤＣＴによるプレシジョンを示す。グラフ９００の横軸は各低次元変換における特徴量（特徴ベクトル）の次元数を示し、縦軸はプレシジョンを示す。なお、図９は、１時間のビデオからランダムに選択された１万フレームの画像と他の１時間のビデオからランダムに選択された１万フレームの画像についてのプレシジョンを示す。この二つの異なる１時間のビデオは、異なった日の同じ時間帯に同じチャンネルで放送された映像であり、オープニング、エンディングの画像等、同一の画像が幾つか含まれている。各画像は３５２×２４０画素の画像であり、閾値θ0及びθ1は０．９に設定されている。
【００８７】
式（５）の上限値ＵＢMSEAは、式（５８）に示すように、式（５９）、（６０）で示される特徴ベクトルの内積で表される。
【数５８】
【数５９】
【数６０】
そこで、ＭＳＥＡを用いた低次元変換方法では、低次元特徴量として式（５９）、（６０）に示す特徴ベクトルを使用し、その内積が閾値θ1以上となる二画像を類似画像として抽出する。
【００８８】
図１０は、ＤＣＴを用いた低次元変換方法を説明するための概略図である。ＤＣＴを用いた低次元変換方法では、まず、各画像の正規化画像を所定のブロックに分割し、二次元ＤＣＴにより周波数成分に変換する。図１０のブロック１０００内の各要素（ａ0、ａ1、ａ2、...）は、二次元ＤＣＴにより変換されたＤＣＴ係数を示す。図１０に示すように、二次元ＤＣＴにより変換されたＤＣＴ係数は、ＭＰＥＧ、ＪＰＥＧ等で用いられるジグザグスキャンにより低周波数成分側から順に並べられる。そして、式（６１）に示すように、低周波数成分側から順に選択されたＤＣＴ係数による特徴ベクトルがＤＣＴによる低次元特徴量となる。
【数６１】
なお、ＤＣＴ係数は正規化画像から生成され、ＤＣ成分は常に０となるため省略される。
【００８９】
ＤＣＴにより変換された周波数成分の全パワーは元の画像のパワーと一致する。つまり、全ＤＣＴ係数を用いて計算したユークリッド距離がＮＣＣ値に対応するので、ＤＣＴ係数の一部（又は全部）を用いて計算したユークリッド距離が（２−２θ0）以下である場合、ＮＣＣ値がθ0以上となることはない。そこで、この方法では、式（６２）に示すように、それぞれのＤＣＴ成分がａi、ｂiである二つの画像の特徴ベクトルのユークリッド距離を求め、ユークリッド距離が（２−２θ1）（θ1≧θ0）以下となる二画像を類似画像として抽出する。
【数６２】
【００９０】
なお、図９に示すグラフ９００では、式（６２）に示すユークリッド距離が（２−２θ1）以下となる画像のペアの集合をＩＩＰＳ'とし、ＮＣＣ値がθ0以上となる画像のペアの集合をＩＩＰＳとしてプレシジョン及びリコールを算出している。
【００９１】
図９に示すように、（特に次元数が低い場合）プレシジョンは必ずしも高くない。プレシジョンは、最も高くても１５％未満であり（次元数が１０２４の場合）、次元数が６４の場合は５〜６％まで落ちる。また、ＤＣＴによるプレシジョンは、本実施形態によるプレシジョンとほとんど同じ値であるが、ＭＳＥＡによるプレシジョンは、極めて低い値となる。
【００９２】
一方、例えば、ＮＣＣの上限値ＵＢについての閾値θ1をＮＣＣ値についての閾値θ0より大きい値にすると（つまり、類似画像を正確にフィルタリングするのではなく、近似的にフィルタリングすると）、リコールは１００％にならなくなるが、プレシジョンを高くすることができる。
【００９３】
図１１（ａ）は、本実施形態及び他の低次元変換方法においてリコールを１００％以下にしたときのリコールとプレシジョンの関係のグラフを示し、図１１（ｂ）は、その拡大図を示す。図１１（ａ）、（ｂ）に示すグラフ１１００、１１１０では、他の低次元変換方法として、ＭＳＥＡと、ＤＣＴと、正規化輝度ヒストグラム（以下、ＮＩＨ（Normalized intensity histograms）と称する）が用いられる。
【００９４】
図１２は、ＮＩＨを用いた低次元変換方法を説明するための概略図である。ＮＩＨを用いた低次元変換方法では、まず、各画像の正規化画像１２００を（２×２、３×３等の）サブ領域１２１０に分割し、各サブ領域について、一定範囲毎に量子化した正規化輝度のヒストグラム１２２０を算出する。そして、式（６３）に示すように、全てのサブ領域のヒストグラムの分布値ａi,j（ｉはサブ領域の番号、ｊは各ヒストグラムの量子化レベルの番号）を順番に並べたベクトルがＮＩＨによる低次元特徴量となる。
【数６３】
この方法では、式（６４）に示すように、各ヒストグラムの分布値がａi,j、ｂi,jである二つの画像の特徴ベクトルのユークリッド距離を求め、ユークリッド距離が（２−２θ1）以下となる二画像を類似画像として抽出する。
【数６４】
【００９５】
また、この方法では、サブ領域の数と正規化輝度の量子化レベルを変更することにより、特徴量の次元数を変更することができる。例えば、サブ領域を２×２にし、正規化輝度の量子化を４レベルにした場合、特徴量の次元数は、２×２×４となる。ＮＩＨを用いて類似画像を抽出する場合、ＮＣＣの上限値は算出されないのでＮＣＣ値が閾値θ0以上になる画像のペアの集合ＩＩＰＳに対するリコールは１００％とならないが、リコール及びプレシジョンが非常に高くなることが知られている。
【００９６】
なお、図１１に示すグラフ１１００では、式（６４）に示すユークリッド距離が（２−２θ1）以下となる画像のペアの集合をＩＩＰＳ'とし、ＮＣＣ値がθ0以上の画像のペアの集合をＩＩＰＳとしてプレシジョン及びリコールを算出している。
【００９７】
図１１（ａ）、（ｂ）に示すグラフ１１００において、「ｐｒｏｐ １６Ｄ」、「ｐｒｏｐ ６４Ｄ」、「ｐｒｏｐ ２５６Ｄ」、「ｐｒｏｐ １０２４Ｄ」のグラフは、本実施形態の低次元変換方法で特徴量次元数がそれぞれ１６、６４、２５６、１０２４のときのリコールとプレシジョンの関係を示す。また、「ＤＣＴ ５Ｄ」、「ＤＣＴ ４４Ｄ」、「ＤＣＴ ２０９Ｄ」のグラフは、ＤＣＴによる低次元変換方法で特徴量次元数がそれぞれ５、４４、２０９のときのリコールとプレシジョンの関係を示す。また、「ＮＩＨ ２×２×４」、「ＮＩＨ ４×４×４」、「ＮＩＨ ４×４×８」のグラフは、ＮＩＨによる低次元変換方法で特徴量次元数がそれぞれ２×２×４、４×４×４、４×４×８のときのリコールとプレシジョンの関係を示す。また、「ＭＳＥＡ １６Ｄ」、「ＭＳＥＡ ６４Ｄ」、「ＭＳＥＡ ５１２Ｄ」のグラフは、ＭＳＥＡによる低次元変換方法で特徴量次元数がそれぞれ１６、６４、５１２のときのリコールとプレシジョンの関係を示す。グラフ１１００、１１１０の横軸はリコールの値を示し、縦軸はプレシジョンの値を示す。
【００９８】
図１１（ａ）、（ｂ）に示すように、本実施形態の低次元変換方法では、特徴量次元数が６４と２５６の場合、リコールを９７％に保ちつつ、プレシジョンが７０％以上となる。特に、リコールが非常に高い（９５％以上）場合、ＤＣＴ、ＮＩＨ、ＭＳＥＡよりプレシジョンが高くなる。
【００９９】
図１３（ａ）は、図１１（ａ）とは異なる画像セットのペア（それぞれ、類似する放送映像から取得した１万画像）についてのリコールとプレシジョンの関係のグラフを示し、図１３（ｂ）は、その拡大図を示す。図１３（ａ）、（ｂ）のグラフ１３００、１３１０では、ほとんど全てのタイプで、図１１（ａ）、（ｂ）のグラフ１１００、１１１０より高いリコールとプレシジョンの組合せが得られる。図１３（ｂ）に示すように、本実施形態の低次元変換方法で特徴量次元数が１０２４又は２５６のときに最も高いパフォーマンスが得られ、次いでＭＳＥＡによる低次元変換方法で特徴量次元数が５１２のとき、本実施形態の低次元変換方法で特徴量次元数が６４又は１６のとき、ＭＳＥＡによる低次元変換方法で特徴量次元数が６４のときの順に高いパフォーマンスが得られる。一方、ＤＣＴによる低次元変換方法では、あまり高いパフォーマンスが得られない。また、ＭＳＥＡによる低次元変換方法で特徴量次元数が１６の場合のパフォーマンスは、本実施形態の低次元変換方法で特徴量次元数が１６の場合より非常に低くなる。つまり、本実施形態の低次元変換方法は、全ての特徴量次元数において、他の低次元変換方法より高いパフォーマンスを得ることができる。
【０１００】
図１４（ａ）、（ｂ）、（ｃ）、（ｄ）は、それぞれリコールを０．９５、０．９７、０．９８、０．９９に固定したときの各低次元変換方法での特徴量次元数とプレシジョンの関係のグラフを示す。図１４（ａ）、（ｂ）、（ｃ）、（ｄ）に示すグラフ１４００、１４１０、１４２０、１４３０において、それぞれ、グラフ１４０１、１４１１、１４２１、１４３１は本実施形態の低次元変換方法でのプレシジョンを示し、グラフ１４０２、１４１２、１４２２、１４３２はＭＳＥＡによる低次元変換方法でのプレシジョンを示し、グラフ１４０３、１４１３、１４２３、１４３３はＤＣＴによる低次元変換方法でのプレシジョンを示し、グラフ１４０４、１４１４、１４２４、１４３４はＮＩＨによる低次元変換方法でのプレシジョンを示す。グラフ１４００〜１４３０の横軸は特徴量次元数を示し、縦軸はプレシジョンを示す。
【０１０１】
図１４（ａ）〜（ｄ）に示すように、リコールを非常に高い値にした場合、本実施形態の低次元変換方法でのプレシジョンは、他の低次元変換方法でのプレシジョンより高くなることが明らかである。図１４（ｄ）に示すように、リコールが０．９９の場合、本実施形態の低次元変換方法では、ＤＣＴによる低次元変換方法よりわずかにパフォーマンスが低くなるが、この場合、そもそもプレシジョンが０．１４以下であり、極めて低い。また、例えば、図１４（ａ）に示すように、リコールが０．９５の場合、本実施形態の低次元変換方法において、特徴量次元数が低い方が高いパフォーマンスを示している（特徴量次元数が６４、１２８、２５６の方が、特徴量次元数が５１２のときよりもプレシジョンが高い）。類似した現象がＭＳＥＡでもみられる（特徴量次元数が１６、３２、６４の方が、特徴量次元数が１２８のときよりもプレシジョンが高い）。これは、これらの方法が複数ブロックへ分割するものであるためと考えられる。そのため、リコールを１００％としない場合、適切なサイズのブロックに分割する必要がある。図１４（ａ）〜（ｄ）から、本実施形態の低次元変換方法でのプレシジョンは、特徴量次元数が６４のとき、すなわち、分割するブロック数が３２のとき、最も高くなると考えられる。
【０１０２】
図１５（ａ）〜（ｄ）は、図１４（ａ）〜（ｄ）とは異なる画像セットのペアについて、それぞれリコールを０．９５、０．９７、０．９８、０．９９に固定したときの各低次元変換方法での特徴量次元数とプレシジョンの関係のグラフを示す。図１５（ａ）、（ｂ）、（ｃ）、（ｄ）に示すグラフ１５００、１５１０、１５２０、１５３０において、それぞれ、グラフ１５０１、１５１１、１５２１、１５３１は本実施形態の低次元変換方法でのプレシジョンを示し、グラフ１５０２、１５１２、１５２２、１５３２はＭＳＥＡによる低次元変換方法でのプレシジョンを示し、グラフ１５０３、１５１３、１５２３、１５３３はＤＣＴによる低次元変換方法でのプレシジョンを示し、グラフ１５０４、１５１４、１５２４、１５３４はＮＩＨによる低次元変換方法でのプレシジョンを示す。グラフ１５００〜１５３０の横軸は特徴量次元数を示し、縦軸はプレシジョンを示す。
【０１０３】
図１５（ａ）〜（ｄ）に示すグラフ１５００〜１５３０では、プレシジョンが比較的高くなっている。例えば、図１５（ｄ）では、リコールが０．９９に固定されているにも関わらず、プレシジョンは略０．６となる。図１５（ａ）〜（ｃ）に示すように、リコールが比較的低い（０．９５〜０．９８）場合、本実施形態の低次元変換方法とＭＳＥＡによる低次元変換方法はＤＣＴ及びＮＩＨによる低次元変換方法よりパフォーマンスがよい。一方、図１５（ｄ）に示すように、リコールが高い（０．９９）場合、本実施形態の低次元変換方法とＤＣＴ及びＮＩＨによる低次元変換方法はＭＳＥＡによる低次元変換方法よりパフォーマンスが高い。つまり、全体として、本実施形態の低次元変換方法は他の低次元変換方法より高いパフォーマンスを得ることができる。
【０１０４】
以上、図９、図１１（ａ）、（ｂ）、図１３（ａ）、（ｂ）、図１４（ａ）〜（ｄ）、図１５（ａ）〜（ｄ）を用いて説明したように、本実施形態の低次元変換方法を用いることにより、他の低次元変換方法より高精度に類似性の高い画像を抽出することが可能となる。
【０１０５】
本実施形態の低次元変換方法では、ラグランジュの未定乗数決定法を用いてＮＣＣの上限値を算出する。これにより、ＮＣＣ値が上限値に近いとき、そのＮＣＣの微分係数は非常に０に近くなり、上限値から離れているとき、微分係数は０から離れる。そのため、ＮＣＣ値が上限値に近い（リコールが高い）場合、上限値をわずかに大きくしても、微分係数が０に近いのでリコールは大きく変化しない。一方、ＮＣＣ値が上限値から離れている場合、微分係数は０から離れているので、プレシジョンは大幅に増加する。これにより、本実施形態の低次元変換方法では、リコールを高く保ちつつ、プレシジョンを高くできると考えられる。
【０１０６】
図１６は、本実施形態による照合処理の時間を示す表である。図１６の表１６００は、それぞれ１万画像（３５２×２４０画素）からなる二つの画像セットにおいて、全ての画像ペアの組合せについて特徴ベクトルのユークリッド距離により全数照合したときの時間を示す。なお、この時間には、特徴ベクトルの算出時間は含まれない。
【０１０７】
表１６００の行１６１０と行１６３０は低次元特徴量による照合処理を、行１６２０と行１６４０は低次元特徴量による照合処理及びそれにより抽出した類似画像のペアに対するＮＣＣによる照合処理を示す。また、行１６１０と行１６２０ではリコールが１００％になるようにし（低次元特徴量による照合処理における閾値θ1＝ＮＣＣによる照合処理における閾値θ0＝０．９）、行１６３０と行１６４０ではリコールが９５％になるようにしている（θ1＞θ0）。各行において、さらに特徴量の次元数毎に処理結果が示される。一方、列１６５０は特徴量次元数を、列１６６０はプレシジョンを、列１６７０はリコールを、列１６８０は処理時間（秒）をそれぞれ示す。
【０１０８】
行１６１０と行１６３０に示すように、照合処理時間は、概ね特徴量の次元数に比例している。これに基づくと、全ての画像ペアの組合せについてＮＣＣによる照合処理（つまり、３５２×２４０次元の演算）を実施した場合、約４０時間かかると考えられる。リコールが１００％の場合、行１６１０に示すようにプレシジョン値は比較的低くなり、類似画像として多くの画像が抽出されるため、行１６２０に示すようにＮＣＣによる照合処理時間は非常に長くなる。一方、リコールが９５％の場合、行１６３０に示すようにプレシジョンが大幅に高くなるため、類似画像の抽出数が抑制され、行１６４０に示すようにＮＣＣによる照合処理時間はあまり長くならない。従って、わずかにリコールを減少させる（９５％）のがリーズナブルであると考えられる。
【０１０９】
なお、上述したように、行１６１０と行１６３０の照合処理時間は、特徴量の次元数に応じて異なるものであり、本実施形態とＭＳＥＡとで差は生じない。また、１万画像（３５２×２４０画素）について、本実施形態に基づく低次元特徴量の算出時間とＭＳＥＡに基づく低次元特徴量の算出時間は、その次元数に関わらず、それぞれ、約３５秒と約３２秒であり、ほとんど差がなかった。つまり、次元数が同じ場合、本実施形態の低次元変換方法とＭＳＥＡによる低次元変換方法とで類似画像を抽出するのにかかる時間は、ほとんど同じとなる。一方、本実施形態の低次元変換方法では、ＭＳＥＡによる低次元変換方法より、類似画像として抽出する画像のペアが少なくなるので、その後のＮＣＣによる照合処理時間は短くなり、トータルの照合時間は短くなる。
【０１１０】
図１７は、画像照合装置を画像符号化装置に適用する場合の例を示す概略構成図である。図１７に示す画像照合装置２は、図１に示す画像照合装置１の各部に加えて画像入力部５３を有する。また、画像照合装置２の制御部６０は、インテグラル画像生成部６１、符号化処理部６２、画像分割部６３、特徴ベクトル算出部６４、上限値算出部６５、第１照合部６６、ＮＣＣ値算出部６７及び第２照合部６８を有する。なお、画像分割部６３は、図１に示す第１画像分割部２２と第２画像分割部２５の機能を備え、特徴ベクトル算出部６４は、図１に示す第１特徴ベクトル算出部２３と第２特徴ベクトル算出部２６の機能を備える。
【０１１１】
画像入力部５３は、ＣＣＤ、ＣＭＯＳ等の光電変換器で構成された２次元検出器と、その２次元検出器上に撮影対象の像を結像する結像光学系等を有する。画像入力部５３は、一定の時間間隔（例えば１／３０秒）毎に撮影を行い、撮影画像を、例えば３５２×２４０画素のデジタル画像に変換し、そのデジタル画像を記憶部５２に記憶する。
【０１１２】
図１８は、図１７に示す画像照合装置２による画像の符号化処理の動作を示すフローチャートである。以下、図１８に示したフローチャートを参照しつつ、画像符号化処理の動作を説明する。なお、以下に説明する動作のフローは、予め記憶部５２に記憶されているプログラムに基づき主に制御部６０により画像照合装置２の各要素と協同して実行される。
【０１１３】
最初に、画像入力部５３は、撮影対象を撮影した画像をデジタル画像に変換し、そのデジタル画像を記憶部５２に記憶する（ステップＳ１８０１）。
【０１１４】
次に、インテグラル画像生成部６１は、記憶部１２に保存されたデジタル画像を読み出し、インテグラル画像を生成し、記憶部５２に記憶する（ステップＳ１８０２）。
【０１１５】
デジタル画像及びインテグラル画像の水平方向をｘ軸、垂直方向をｙ軸とし、デジタル画像の座標（ａ、ｂ）における正規化画素値をＩ（ａ、ｂ）とする。インテグラル画像生成部６７は、各画素が式（６５）、（６６）からなるインテグラル画像（Viola, P., Jones, M., May 2004. Robust real-time face detection. International Journal of Computer Vision 57 (2), 137-154.）をそれぞれ生成する。
【数６５】
【数６６】
つまり、式（６５）により算出されるインテグラル画像の各画素値は、原点からその画素までの全ての画素の正規化画素値の総和であり、式（６６）により算出されるインテグラル画像の各画素値は、原点からその画素までの全ての画素の正規化画素値の二乗和である。
【０１１６】
図１９は、インテグラル画像を説明するための概略図である。例えば、インテグラル画像１９００において、座標（ｘ1、ｙ1）、（ｘ1、ｙ2）、（ｘ2、ｙ1）、（ｘ2、ｙ2）で囲まれた領域Ｄの正規化画素値の総和と正規化画素値の二乗和は、それぞれ、式（６７）、（６８）により算出することができる。
【数６７】
【数６８】
従って、予めインテグラル画像を作成しておくことにより、所定領域内の正規化画素値の総和と二乗和の算出処理を高速化できるので、特徴ベクトルの算出を高速化できる。
【０１１７】
次に、符号化処理部６２は、記憶部１２に保存されたデジタル画像を読み出し、デジタル画像のフォーマット（画素数）変換、８×８画素、１６×１６画素等の符号化ブロックへの分割等の符号化前処理を実施する（ステップＳ１８０３）。
【０１１８】
以下のステップＳ１８０４〜Ｓ１８１４の処理は、符号化ブロック毎に実施される。まず、画像分割部６３は、符号化ブロックをさらに所定サイズの複数のブロックに分割する（ステップＳ１８０４）。例えば、画像分割部６３は、１６×１６画素の画像を４×４画素の１６ブロックに分割する。
【０１１９】
次に、特徴ベクトル算出部６４は、インテグラル画像を用いて、符号化ブロックについて、式（３７）に示す特徴ベクトルを算出し、符号化ブロックと関連付けて記憶部１２に保存する（ステップＳ１８０５）。
【０１２０】
次に、上限値算出部６５は、記憶部１２から照合元ブロックの特徴ベクトルを取得する（ステップＳ１８０６）。
【０１２１】
この照合元ブロックは、例えば、１フレーム前のデジタル画像における各符号化ブロックである。なお、前方向予測だけでなく、後方向予測も用いる場合、上限値算出部６５は、１フレーム後のデジタル画像における各符号化ブロックも照合元ブロックとして用いる。
【０１２２】
次に、上限値算出部６５は、符号化ブロックの特徴ベクトルと照合元ブロックの特徴ベクトルの内積を算出して、符号化ブロックと照合元ブロックについて、ＮＣＣの上限値ＵＢを求める（ステップＳ１８０７）。
【０１２３】
次に、第１照合部６６は、符号化ブロックと照合元ブロックについての上限値ＵＢが閾値θ1以上であるか否かを判定する（ステップＳ１８０８）。
【０１２４】
上限値ＵＢが閾値θ1以上である場合、第１照合部６５は、その照合元ブロックを符号化ブロックと類似すると判定し、ＮＣＣ値算出部６７は、その符号化ブロックと照合元ブロックについてＮＣＣ値を算出する（ステップＳ１８０９）。
【０１２５】
次に、第２照合部６８は、ＮＣＣ値算出部６７が算出したＮＣＣ値がその符号化ブロックについて算出されたＮＣＣ値のうち最大であるか否かを判定する（ステップＳ１８１０）。
【０１２６】
ＮＣＣ値算出部６７が算出したＮＣＣ値がその符号化ブロックについて算出されたＮＣＣ値のうち最大である場合、第２照合部６８は、その照合元ブロックを動き検出用ブロックとして記憶部５２に記憶（すでに記憶されている場合、更新）する（ステップＳ１８１１）。
【０１２７】
一方、上限値ＵＢが閾値θ1未満である場合、ＮＣＣ値算出部６７が算出したＮＣＣ値がその符号化ブロックについて算出されたＮＣＣ値のうち最大でない場合、又は、第２照合部６８が照合元ブロックを動き検出用ブロックとして記憶部５２に記憶した場合、制御部６０は、符号化ブロックを全ての照合元ブロックと比較したか否かを判定する（ステップＳ１８１２）。
【０１２８】
制御部６０は、まだ比較していない照合元ブロックがある場合、ステップＳ１８０７〜Ｓ１８１１の処理を繰り返し、全ての照合元ブロックとの比較が完了すると、符号化処理部６２は、記憶部５２に記憶された動き検出用ブロックの候補を用いて動き補償を行う。また、符号化処理部６２は、ＤＣＴ変換、量子化等の符号化処理を実施する（ステップＳ１８１３）。
【０１２９】
次に、制御部６０は、全ての符号化ブロックの符号化処理が完了したか否かを判定する（ステップＳ１８１４）。全ての符号化ブロックの符号化処理が完了していない場合、ステップＳ１８０３〜Ｓ１８１３の処理を繰り返し、全ての符号化ブロックの符号化処理が完了すると、制御部６０は、一連のステップを終了する。
【０１３０】
このようにして符号化されたデータは、インターフェース部５１を介して、外部のデコード装置（不図示）に送信される。
【０１３１】
なお、第２照合部６８は、符号化ブロックと全ての照合元ブロックについての上限値ＵＢが閾値θ1未満であった場合は、その符号化ブロックに類似する照合元ブロックが存在しないと判断し、任意の照合元ブロックを動き検出用ブロックとして選択する。あるいは、その場合、上限値ＵＢが最も高かった照合元ブロックを動き検出用ブロックとして選択してもよい。
【０１３２】
以上詳述したように、図１８に示したフローチャートに従って動作することによって、画像照合装置２は、画像符号化装置として動作し、高速に動き検出処理を実施することが可能となった。また、インテグラル画像を用いて特徴量算出処理を高速化することが可能となった。また、本実施形態による方法は、スライディングウィンドウに基づく探索（Wei, S.-D., Lai, S.-H., 2007. Efficient normalized cross correlation based on adaptive multilevel successive elimination. In: Proceedings of the 8th Asian conference on Computer vision - Volume Part I. ACCV’07. Springer-Verlag, Berlin, Heidelberg, pp. 638-646.）とも親和性が高い。
【符号の説明】
【０１３３】
１、２ 画像照合装置
１１、５１ インターフェース部
１２、５２ 記憶部
２０、３０、４０、６０ 制御部
２１ 第１画像取得部
２２ 第１画像分割部
２３ 第１特徴ベクトル算出部
２４ 第２画像取得部
２５ 第２画像分割部
２６ 第２特徴ベクトル算出部
２７、６５ 上限値算出部
２８、６６ 第１照合部
３１、６７ ＮＣＣ値算出部
３２、６８ 第２照合部
４１ インデックス判定部
５３ 画像入力部
６１ インテグラル画像生成部
６２ 符号化処理部
６３ 画像分割部
６４ 特徴ベクトル算出部
【技術分野】
【０００１】
本発明は、画像照合装置、画像照合方法及びコンピュータプログラムに関し、特に、二つの画像が類似するか否かを判定する画像照合装置、画像照合方法及びコンピュータプログラムに関する。
【背景技術】
【０００２】
画像照合技術は、テンプレートマッチング、ブロック動き補償、画像圧縮、ステレオビジョン等の、多くの画像を用いたアプリケーションにおける重要な技術の一つである。最近は、画像の近似コピー検出、画像マイニング及びその画像アノテーションへの応用等、大規模な画像／映像データベースにおける画像照合技術の新しいアプリケーションについても研究が開始されている。
【０００３】
一般に、画像照合の方法として、画素値の差の絶対値の合計（ＳＡＤ：Sum of Absolute Difference）、画素値の差の二乗の合計（ＳＳＤ：Sum of Squared Difference）、正規化相互相関（ＮＣＣ：Normalized Cross-Correlation）等が用いられている。特に、ＮＣＣは、画像間の輝度のずれ、コントラストの違い等による影響が少なく、輝度の補正等がなされた画像間の照合にも適している。例えば、水平方向をｘ軸、垂直方向をｙ軸とする画像Ｉ1と画像Ｉ2についてのＮＣＣ値は、式（１）で定義される。
【数１】
ここで、φx,yは画像Ｉ1の座標（ｘ、ｙ）の画素の画素値であり、ψx,yは画像Ｉ2の座標（ｘ、ｙ）の画素の画素値であり、φaveは画像Ｉ1の全画素の画素値の平均値であり、ψaveは画像Ｉ2の全画素の画素値の平均値である。式（１）に示すように、ＮＣＣ値を算出するためには、対応する画素毎に画素値の相関を計算する必要があるため、特に、照合する画像のサイズが大きい場合には、計算量が大きくなり、算出処理に多大な時間が必要となる。
【０００４】
そこで、例えば、非特許文献１には、ＮＣＣに基づくテンプレートマッチング処理を高速化するための技術が開示されている。非特許文献１に開示された技術は、画像全体に対してテンプレートをずらしながらスキャンする場合、畳み込み演算を用いて各位置での照合結果を求める。一方、空間領域間の畳み込み演算は、周波数領域間での単純な掛け算に変換できるので、空間領域のデータを高速フーリエ変換（ＦＦＴ：Fast Fourie Transform）を用いて周波数領域に変換することにより処理を高速化できる。しかし、この技術はある画像に対してテンプレートをスキャンする場合には効果的であるが、複数の画像に対して照合をする場合は各画像に対してＦＦＴ及び逆ＦＦＴを行う必要があり、かえって処理が遅くなるおそれがある。
【０００５】
一方、非特許文献２には、ＮＣＣの演算処理を効率化するための技術として、ＭＳＥＡ（multilevel successive elimination algorithm）が開示されている。非特許文献２に開示された技術では、式（２）に示すCauchy-Schwarzの不等式に基づいてＮＣＣの上限値が算出される。
【数２】
【０００６】
画像Ｉ1と画像Ｉ2についてブロック内画素数がｍである、ｎ個のブロックに分割した場合、式（１）のＮＣＣ値は、ブロックの番号ｉ（１≦ｉ≦ｎ）及びブロック内の画素の番号ｊ（１≦ｊ≦ｍ）を用いて式（３）で表される。
【数３】
ここで、ｘi,j及びｙi,jは各画素の正規化画素値であり、式（４）で定義される。
【数４】
式（２）、（３）から、ＮＣＣの上限値ＵＢMSEAは、式（５）で表される。
【数５】
ここで、ＸＸi及びＹＹiは式（６）で定義される。
【数６】
ＮＣＣ値はこの上限値ＵＢMSEAを越えることがないので、ＮＣＣ値の代わりに、ＮＣＣ値よりも簡易に算出できる上限値ＵＢMSEAを用いて二つの画像の類似性を評価することができる。
【先行技術文献】
【非特許文献】
【０００７】
【非特許文献１】Tsai, D.-M., Lin, C.-T., 2003. Fast normalized cross correlation for defect detection. Pattern Recognition Letters 24 (15), 2625-2631.
【非特許文献２】Wei, S.-D., Lai, S.-H., 2007. Efficient normalized cross correlation based on adaptive multilevel successive elimination. In: Proceedings of the 8th Asian conference on Computer vision - Volume Part I. ACCV’07. Springer-Verlag, Berlin, Heidelberg, pp. 638-646
【発明の概要】
【発明が解決しようとする課題】
【０００８】
非特許文献２に開示された技術を用いることにより、高速に画像照合を実施することができる。しかしながら、非特許文献２に開示されたＮＣＣの上限値は、ＮＣＣ値より大きくなりすぎる場合がある。その場合、照合する二画像のＮＣＣ値が低いにも関わらず、その二画像の類似性が高いと判定するおそれがあった。
【０００９】
そこで、本発明の目的は、画像照合を高速に実施するとともに、画像の類似性を高精度に判定することが可能な画像照合装置、画像照合方法及びそのような画像照合方法をコンピュータに実行させるコンピュータプログラムを提供することにある。
【課題を解決するための手段】
【００１０】
本発明に係る画像照合装置は、第１画像を複数のブロックに分割し、第２画像を第１画像と同数のブロックに分割する画像分割部と、第１画像及び第２画像のそれぞれについて、各ブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づき、第１画像の特徴ベクトル及び第２画像の特徴ベクトルを算出する特徴ベクトル算出部と、ラグランジュの未定乗算決定法を利用して２つの画像の正規化相互相関値を求める式から導いた２つの画像の正規化相互相関の上限値を算出する式を用いて、算出された第１画像の特徴ベクトル及び第２画像の特徴ベクトルに基づき、第１画像と第２画像の正規化相互相関の上限値を算出する上限値算出部と、上限値が第１の閾値以上であるか否かに基づいて第１画像と第２画像を照合する第１照合部と、を有する。
【００１１】
さらに、本発明に係る画像照合装置において、特徴ベクトル算出部は、ｍをブロック内の画素数、ｎをブロックの数、φi,jを第１画像におけるｉ番目のブロック内のｊ番目の画素の画素値、φaveを第１画像の全画素の画素値の平均値として、次のベクトルξxを第１画像の特徴ベクトルとして算出し、
【数７】
ψi,jを第２画像におけるｉ番目のブロック内のｊ番目の画素の画素値、ψaveを第２画像の全ての画素の画素値の平均値として、次のベクトルξyを第２画像の特徴ベクトルとして算出し、
【数８】
上限値算出部は、第１画像の特徴ベクトルと第２画像の特徴ベクトルの内積を上限値として算出することが好ましい。
【００１２】
さらに、本発明に係る画像照合装置において、上限値が第１の閾値以上である第１画像と第２画像の正規化相互相関値を算出する正規化相互相関値算出部と、正規化相互相関値が第１の閾値以下の値である第２の閾値以上であるか否かに基づいて第１画像と第２画像を詳細に照合する第２照合部と、をさらに有することが好ましい。
【００１３】
また、本発明に係る画像照合方法は、第１画像を複数のブロックに分割するステップと、第１画像について、各ブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づき、第１画像の特徴ベクトルを算出するステップと、第２画像を第１画像と同数のブロックに分割するステップと、第２画像について、各ブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づき、第２画像の特徴ベクトルを算出するステップと、ラグランジュの未定乗算決定法を利用して２つの画像の正規化相互相関値を求める式から導いた２つの画像の正規化相互相関の上限値を算出する式を用いて、算出された第１画像の特徴ベクトル及び第２画像の特徴ベクトルに基づき、第１画像と第２画像の正規化相互相関の上限値を算出するステップと、上限値が第１の閾値以上であるか否かに基づいて第１画像と第２画像を照合するステップと、を含む。
【００１４】
さらに、本発明に係る画像照合方法において、第１画像の特徴ベクトルを算出するステップにおいて、ｍをブロック内の画素数、ｎをブロックの数、φi,jを第１画像におけるｉ番目のブロック内のｊ番目の画素の画素値、φaveを第１画像の全画素の画素値の平均値として、次のベクトルξxを第１画像の特徴ベクトルとして算出し、
【数９】
第２画像の特徴ベクトルを算出するステップにおいて、ｍをブロック内の画素数、ｎをブロックの数、ψi,jを第２画像におけるｉ番目のブロック内のｊ番目の画素の画素値、ψaveを第２画像の全ての画素の画素値の平均値として、次のベクトルξyを第２画像の特徴ベクトルとして算出し、
【数１０】
上限値を算出するステップにおいて、第１画像の特徴ベクトルと第２画像の特徴ベクトルの内積を上限値として算出することが好ましい。
【００１５】
さらに、本発明に係る画像照合方法において、上限値が第１の閾値以上である第１画像と第２画像の正規化相互相関値を算出するステップと、正規化相互相関値が第１の閾値以下の値である第２の閾値以上であるか否かに基づいて第１画像と第２画像を詳細に照合するステップと、をさらに有することが好ましい。
【００１６】
また、本発明に係るコンピュータプログラムは、第１画像を複数のブロックに分割するステップと、第１画像について、各ブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づき、第１画像の特徴ベクトルを算出するステップと、第２画像を第１画像と同数のブロックに分割するステップと、第２画像について、各ブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づき、第２画像の特徴ベクトルを算出するステップと、ラグランジュの未定乗算決定法を利用して２つの画像の正規化相互相関値を求める式から導いた２つの画像の正規化相互相関の上限値を算出する式を用いて、算出された第１画像の特徴ベクトル及び第２画像の特徴ベクトルに基づき、第１画像と第２画像の正規化相互相関の上限値を算出するステップと、上限値が第１の閾値以上であるか否かに基づいて第１画像と第２画像を照合するステップと、をコンピュータに実行させる。
【発明の効果】
【００１７】
本発明によれば、画像照合を高速に実施するとともに、画像の類似性を高精度に判定することが可能な画像照合装置、画像照合方法及びそのような画像照合方法をコンピュータに実行させるコンピュータプログラムを提供することができる。
【図面の簡単な説明】
【００１８】
【図１】本発明を適用した画像照合装置の概略構成図である。
【図２】画像照合処理の動作を示すフローチャートである。
【図３】ブロック毎に分割した画像の例である。
【図４】制御部の他の例を示す概略構成図である。
【図５】詳細な画像照合処理の動作を示すフローチャートである。
【図６】制御部の他の例を示す概略構成図である。
【図７】画像の取得処理の動作を示すフローチャートである。
【図８】画像の照合処理の動作を示すフローチャートである。
【図９】特徴量次元数毎のプレシジョンを表すグラフである。
【図１０】ＤＣＴを用いた低次元特徴量を説明するための概略図である。
【図１１】（ａ）、（ｂ）は、リコールとプレシジョンの関係のグラフである。
【図１２】ＮＩＨを用いた低次元特徴量を説明するための概略図である。
【図１３】（ａ）、（ｂ）は、リコールとプレシジョンの関係のグラフである。
【図１４】（ａ）〜（ｄ）は、特徴量次元数毎のプレシジョンを表すグラフである。
【図１５】（ａ）〜（ｄ）は、特徴量次元数毎のプレシジョンを表すグラフである。
【図１６】特徴量次元数毎の照合処理時間を示す表である。
【図１７】画像照合装置の他の例を示す概略構成図である。
【図１８】画像の符号化処理の動作を示すフローチャートである。
【図１９】インテグラル画像を説明するための概略図である。
【発明を実施するための形態】
【００１９】
以下、本発明に係る画像照合装置、画像照合方法及びコンピュータプログラムについて図を参照しつつ説明する。但し、本発明の技術的範囲はそれらの実施の形態に限定されず、特許請求の範囲に記載された発明とその均等物に及ぶ点に留意されたい。
【００２０】
図１は、本発明を適用した画像照合装置の概略構成を示す図である。図１に示すように、画像照合装置１は、インターフェース部１１、記憶部１２、表示部１３及び制御部２０を有する。以下、画像照合装置１の各部について詳細に説明する。
【００２１】
インターフェース部１１は、例えばインターネット、電話回線網（携帯端末回線網、一般電話回線網を含む）、イントラネット等のネットワークを介して他のコンピュータ等と画像データ及び各種のデータを送受信する通信インターフェースであり、接続するネットワークの通信インターフェース回路を有する。また、インターフェース部１１は、例えばＵＳＢ等のシリアルバスに準じるインターフェース回路を有し、フラッシュメモリ等を接続し、そのフラッシュメモリ等から画像データ及び各種のデータを取得するようにしてもよい。インターフェース部１１は、制御部２０と接続されており、制御部２０により制御される。
【００２２】
記憶部１２は、ＲＡＭ、ＲＯＭ等のメモリ装置、ハードディスク等の固定ディスク装置、又はフレキシブルディスク、光ディスク等の可搬用の記憶装置等を有する。また、記憶部１２には、画像照合装置１の各種処理に用いられるコンピュータプログラム、データベース、テーブル等が格納される。記憶部１２は、制御部２０と接続され、インターフェース部１１を介して取得した画像データを格納するとともに、制御部２０により画像データについてなされた演算結果を格納する。
【００２３】
制御部２０は、複数の画像について特徴ベクトルを算出し、各画像が類似するか否かを判定する。そのために、制御部２０は、第１画像取得部２１、第１画像分割部２２、第１特徴ベクトル算出部２３、第２画像取得部２４、第２画像分割部２５、第２特徴ベクトル算出部２６、上限値算出部２７及び第１照合部２８を有する。また、制御部２０は、インターフェース部１１、記憶部１２及び表示部１３と接続され、インターフェース部１１のデータ送受信制御、記憶部１２の制御、表示部１３の表示制御等を行う。制御部２０は、予め記憶部１２に記憶されているプログラムに基づいて動作する。あるいは、制御部２０は、集積回路、マイクロプロセッサ、ファームウェア等で構成されてもよい。
【００２４】
図２は、画像照合装置１による画像照合処理の動作を示すフローチャートである。以下、図２に示したフローチャートを参照しつつ、画像照合処理の動作を説明する。なお、以下に説明する動作のフローは、予め記憶部１２に記憶されているプログラムに基づき主に制御部２０により画像照合装置１の各要素と協同して実行される。
【００２５】
最初に、第１画像取得部２１は、インターフェース部１１を介して、外部のコンピュータ、フラッシュメモリ等から画像を取得し（以下、第１画像取得部２１が取得した画像を第１画像と称する）、記憶部１２に保存する（ステップＳ２０１）。
【００２６】
次に、第１画像分割部２２は、記憶部１２に保存された第１画像を読み出し、第１画像を所定サイズの複数のブロックに分割する（ステップＳ２０２）。例えば、第１画像分割部２２は、３５２×２４０画素の画像を２２×１５画素の２５６ブロックに分割する。
【００２７】
図３は、ブロック毎に分割した画像の例を示す模式図である。図３の例では、画像３００を、各ブロックの画素数がｍであるｎ個のブロックに分割した画像３１０を示す。
【００２８】
次に、第１特徴ベクトル算出部２３は、第１画像について、分割したブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づく特徴ベクトルを算出し、第１画像と関連付けて記憶部１２に保存する（ステップＳ２０３）。
【００２９】
画像ｘ、画像ｙについてブロック内画素数がｍである、ｎ個のブロックに分割した場合のＮＣＣ値を、ブロックの番号ｉ（１≦ｉ≦ｎ）及びブロック内の画素の番号ｊ（１≦ｊ≦ｍ）を用いて式（１１）で定義する。
【数１１】
ここで、φi,j、ψi,jは、各画素の画素値であり、φaveは画像ｘの全画素の画素値の平均値であり、ψaveは画像ｙの全画素の画素値の平均値である。
【００３０】
画像ｘ及び画像ｙの各画素の正規化画素値ｘi,j及びｙi,jを式（１２）で定義する。
【数１２】
なお、正規化画素値は、厳密には式（１３）とすべきであるが、本実施形態では単純化のため正規化画素値を式（１２）で定義する。
【数１３】
【００３１】
式（１１）のＮＣＣ値は、式（１２）の正規化画素値ｘi,j及びｙi,jを用いて式（１４）で表される。
【数１４】
【００３２】
また、画像ｘのｉ番目のグループにおける正規化画素値ｘi,jの総和及び二乗和並びに画像ｙのｉ番目のグループにおける正規化画素値ｙi,jの総和及び二乗和を式（１５）〜（１８）で定義する。
【数１５】
【数１６】
【数１７】
【数１８】
ここで、正規化画素値の性質上、式（１９）、（２０）が成り立つ。
【数１９】
【数２０】
【００３３】
以下、ラグランジュの未定乗数決定法を利用して式（１４）から２つの画像のＮＣＣの厳密な上限値を算出する式を求める。本実施形態では、式（２１）の関数の極値問題を考える。
【数２１】
ここで、λ1,i、λ2,i、λ3,i、λ4,iは、４ｎ個のラグランジュ未定乗数である。式（２１）において、微分係数ｄΛ＝０（全てのｉ、ｊについて∂Λ/∂ｘi,j＝０、∂Λ/∂yi,j＝０）とすることにより、式（２２）、（２３）の２ｍｎ個の式が得られる。
【数２２】
【数２３】
式（２２）、（２３）から式（２４）が得られる。
【数２４】
【００３４】
ここで、式（２４）が成立すると仮定すると、ｘi,j、ｙi,jは式（２５）で求められる。
【数２５】
この場合、式（２５）の右辺はｉ（ブロックの番号）のみにより定まるので、ｘi,j、ｙi,jは全てのｊ（ブロック内の画素）について同一の値をとることになり、条件が強すぎることになる。一方、式（２４）が成立しないと仮定すると、式（２４）の左辺の行列の行列式から式（２６）が得られる。
【数２６】
また、式（２２）、（２３）、（２６）から、式（２７）が得られる。
【数２７】
この場合、式（２２）と式（２３）は、同じ意味を表すことになる。従って、本実施形態では、式（２２）のみを用いてＮＣＣの上限値を算出する式を求める。
【００３５】
まず、式（２２）から式（２８）が得られる。
【数２８】
式（２８）を式（１７）に代入することにより、式（２９）が得られる。
【数２９】
一方、式（２８）を式（１８）に代入することにより、式（３０）が得られる。
【数３０】
式（２９）より、式（３１）が得られる。
【数３１】
式（３１）を式（３０）に代入することにより、式（３２）が得られる。
【数３２】
式（２８）、（３１）、（３２）を式（１４）に代入することにより、式（３３）が得られる。
【数３３】
【００３６】
つまり、ＮＣＣの上限値ＵＢは式（３４）で与えられ、ＮＣＣの下限値ＬＢは式（３５）で与えられる。
【数３４】
【数３５】
ＮＣＣ値はこの上限値ＵＢを越えることがないので、ＮＣＣ値の代わりに上限値ＵＢを用いて画像ｘと画像ｙの類似性を評価することができる。
【００３７】
一方、式（３４）の上限値ＵＢは、式（３６）に示されるように、式（３７）、（３８）で示される特徴ベクトルξxとξyの内積となる。
【数３６】
【数３７】
【数３８】
式（３７）、（３８）に示されるように、特徴ベクトルξx及びξyは、それぞれ画像ｘ、画像ｙのみから算出することができる。つまり、画像照合装置１は、予め画像ｘ、画像ｙをそれぞれ２ｎ個の要素をもつ特徴ベクトルξx及びξyに変換しておくことにより、その内積を算出して上限値ＵＢを求めることができ、特に、ある画像を複数の画像と照合するときに高速に照合することができる。そこで、第１特徴ベクトル算出部２３は、第１画像の特徴ベクトルとして第１画像について式（３７）の特徴ベクトルξxを算出する。
【００３８】
次に、第２画像取得部２４は、インターフェース部１１を介して、外部のコンピュータ、フラッシュメモリ等から画像を取得し（以下、第２画像取得部２４が取得した画像を第２画像と称する）、記憶部１２に保存する（ステップＳ２０４）。
【００３９】
次に、第２画像分割部２５は、記憶部１２に保存された第２画像を読み出し、第１画像を分割したブロックとそれぞれ同サイズ、同数の複数のブロックに分割する（ステップＳ２０５）。
【００４０】
次に、第２特徴ベクトル算出部２６は、第２画像について、分割したブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づく特徴ベクトル、つまり式（３８）の特徴ベクトルξyを算出し、第２画像と関連付けて記憶部１２に保存する（ステップＳ２０６）。
【００４１】
次に、上限値算出部２７は、式（３６）に示すように特徴ベクトルξx、ξyの内積を算出して、上限値ＵＢを求める（ステップＳ２０７）。
【００４２】
次に、第１照合部２８は、上限値ＵＢが所定の閾値θ1以上であるか否かに基づいて第１画像と前記第２画像を照合する。まず、第１照合部２８は、上限値ＵＢが閾値θ1以上であるか否かを判定する（ステップＳ２０８）。この閾値θ1は、ＮＣＣを用いて二つの画像が類似するか否かを判定するときに、ＮＣＣ値と比較する閾値θ0（−１≪θ0＜１）と同じ値である。あるいは、閾値θ1は、閾値θ0より大きい値としてもよい。
【００４３】
第１照合部２８は、上限値ＵＢが閾値θ1以上である場合、第１画像が第２画像と類似すると判定し（ステップＳ２０９）、一連のステップを終了する。一方、第１照合部２８は、上限値ＵＢが閾値θ1未満である場合、第１画像が第２画像と類似しないと判定し（ステップＳ２１０）、一連のステップを終了する。
【００４４】
例えば、画像照合装置１は、類似すると判定した画像のペアを利用者が確認できるように、表示部１３に表示する。あるいは、画像照合装置１は、類似すると判定した画像のペアの情報をインターフェース部１１を介して外部のコンピュータ（不図示）に通知してもよい。
【００４５】
なお、上限値ＵＢと下限値ＬＢの間の範囲が狭いほど、上限値ＵＢは、ＮＣＣ値をより正確に近似していることになる。上限値ＵＢと下限値ＬＢの間の範囲は、式（３４）、（３５）から、式（３９）により求めることができる。
【数３９】
一方、
【数４０】
であり、ここで、Ｇj,k＝ｘi,jｘi,j−ｘi,jｘi,kとすると、式（４１）が成り立つ。
【数４１】
従って、式（４０）は、式（４２）で表される。
【数４２】
ここで、Ｅkはｋにわたる期待値であり、Ｖａｒi(ｘi,j)はｉ番目のブロック内のｘi,jの分散である。同様に、
【数４３】
が成立し、式（３９）は、式（４４）で表される。
【数４４】
【００４６】
つまり、上限値ＵＢと下限値ＬＢの間の範囲を狭くするためには、各ブロック内部の正規化画素値の分散を小さくするか、一ブロック内の画素数ｍを小さくする必要がある。しかし、一ブロック内の画素数ｍを小さくすると、ブロック数ｎが大きくなり、特徴ベクトルξxとξyの次元数が大きくなるため、現実的でない。従って、同一ブロック内の画素の間の正規化画素値の差が可能な限り小さくなるように各ブロックを構成し、各ブロック内部の正規化画素値の分散を小さくすることが好ましい。一般に、画像のマルコフ性により、近傍の画素の（正規化）画素値は、非常に近い値をもつことが知られている。従って、例えばブロックの形状を、長方形、特に正方形に近い形にすることにより、各ブロック内の画素の画素値を近付け、上限値ＵＢをよりＮＣＣ値に近似させることができる。
【００４７】
なお、ＭＳＥＡにより算出されるＮＣＣの上限値ＵＢMSEAは、式（５）から式（４５）のように表される。
【数４５】
本実施形態により算出されるＮＣＣの上限値ＵＢと、ＭＳＥＡにより算出されるＮＣＣの上限値ＵＢMSEAについて、式（３４）と式（４５）のルート内の式を比較することにより式（４６）が成立する。
【数４６】
従って、式（４７）が成立し、本実施形態により算出されるＮＣＣの上限値ＵＢは、ＭＳＥＡにより算出されるＮＣＣの上限値ＵＢMSEA以下となる。
【数４７】
つまり、本実施形態により算出されるＮＣＣの上限値ＵＢは、ＭＳＥＡにより算出されるＮＣＣの上限値ＵＢMSEAより正確にＮＣＣ値を近似していることが理論的に示される。
【００４８】
以上詳述したように、図２に示したフローチャートに従って動作することによって、画像照合装置１は、元の情報量（数万程度の画素数）から大幅に圧縮した（数十〜百程度の）特徴ベクトル（すなわち低次元特徴量）を用いて、画像の類似性を高速かつ高精度に判定することが可能となった。
【００４９】
なお、画像照合装置１は、照合する二画像のサイズが異なる場合、二画像のサイズを同一サイズに正規化することにより、サイズの異なる二画像についても低次元特徴量を用いた画像照合をすることができる。
【００５０】
図４は、制御部の他の例を示す概略構成図である。図４に示す制御部３０は、図１に示す制御部２０の代わりに用いることが可能である。図４に示す制御部３０は、図１に示す制御部２０の各部に加えて、ＮＣＣ値算出部３１と第２照合部３２を有する。
【００５１】
図５は、図４に示す制御部３０を用いる画像照合装置１による詳細な画像照合処理の動作を示すフローチャートである。以下、図５に示したフローチャートを参照しつつ、詳細な画像照合処理の動作を説明する。なお、以下に説明する動作のフローは、予め記憶部１２に記憶されているプログラムに基づき主に制御部３０により画像照合装置１の各要素と協同して実行される。また、図５に示したフローチャートは、図２に示したフローチャートにより、類似すると判定された画像のペアに対して実行される。
【００５２】
最初に、ＮＣＣ値算出部３１は、類似すると判定された第１画像と第２画像について式（１１）を用いてＮＣＣ値を算出する（ステップＳ５０１）。
【００５３】
次に、第２照合部３２は、ＮＣＣ値算出部３１により算出されたＮＣＣ値が所定の閾値θ0（例えば０．９）以上であるか否かに基づいて第１画像と第２画像を詳細に照合する。まず、第２照合部３２は、ＮＣＣ値が閾値θ0以上であるか否かを判定する（ステップＳ５０２）。なお、上述した通り、この閾値θ0は閾値θ1以下の値の閾値である。
【００５４】
第２照合部３２は、ＮＣＣ値が閾値θ0以上である場合、第１画像と第２画像は類似すると判定し（ステップＳ５０３）、一連のステップを終了する。一方、第２照合部２３は、ＮＣＣ値が閾値θ0未満である場合、第１画像と第２画像は類似しないと判定し（ステップＳ５０４）、一連のステップを終了する。
【００５５】
例えば、画像照合装置１は、図５に示したフローチャートにより類似すると判定した画像のペアを表示部１３に表示する。あるいは、画像照合装置１は、図５に示したフローチャートにより類似すると判定した画像のペアの情報をインターフェース部１１を介して外部のコンピュータ（不図示）に通知してもよい。
【００５６】
なお、制御部２０は、第１画像取得部２１と第２画像取得部２４の代わりに、第１画像及び第２画像の何れをも取得できる画像取得部を備えることもできる。また、第１画像分割部２２と第２画像分割部２５の代わりに第１画像及び第２画像の何れをも所定のブロックに分割できる画像分割部を備えることもできる。さらに、第１特徴ベクトル算出部２３と第２特徴ベクトル算出部２６の代わりに第１画像の特徴ベクトル及び第２画像の特徴ベクトルの何れをも算出できる特徴ベクトル算出部を備えることもできる。
【００５７】
また、ステップＳ２０４〜Ｓ２０６の第２画像に対する処理は、ステップＳ２０１〜Ｓ２０３の第１画像に対する処理より前に実施してもよい。また、画像照合装置１が複数のＣＰＵを備えること等により、並列処理が可能な場合には、ステップＳ２０４〜Ｓ２０６の第２画像に対する処理と、ステップＳ２０１〜Ｓ２０３の第１画像に対する処理を並列に実施してもよい。
【００５８】
以上詳述したように、図５に示したフローチャートに従って動作することによって、画像照合装置１は、低次元空間での類似性が高い（ＮＣＣの上限値ＵＢが閾値θ1以上の）画像のペアのうち、オリジナル空間での類似性が低い（ＮＣＣ値が閾値θ0未満の）画像のペアを取り除く（クレンジングする）ことができる。これにより、画像照合装置１は、ある画像のペアに対して低負荷で類似性を判定しておき、その判定において類似性が高いと判定した画像のペアのみについて高負荷で高精度に類似性を判定するので、効率よく画像を照合できるようになった。
【００５９】
図６は、画像照合装置を画像検索サーバに適用する場合の制御部の概略構成図である。図６に示す制御部４０は、図１に示す制御部２０又は図４に示す制御部３０の代わりに用いることが可能である。図６に示す制御部４０は、図４に示す制御部３０の各部に加えて、インデックス判定部４１を有する。なお、以下では、第１画像を事前に記憶部１２に格納された照合元の画像とし、第２画像を外部のコンピュータ（不図示）等から検索を要求された問合せ画像として説明する。
【００６０】
図７は、図６に示す制御部４０を用いる画像照合装置１による画像の取得処理の動作を示すフローチャートである。以下、図７に示したフローチャートを参照しつつ、画像取得処理の動作を説明する。なお、以下に説明する動作のフローは、予め記憶部１２に記憶されているプログラムに基づき主に制御部４０により画像照合装置１の各要素と協同して実行される。
【００６１】
図７に示すフローチャートでは、画像照合装置１は、第１画像の取得のみを実施し、第１画像と第２画像の照合処理は、後述するフローチャートで説明する。図７に示すステップＳ７０１〜Ｓ７０３の処理は、図２に示すステップＳ２０１〜Ｓ２０３の処理と同じであるため、説明を省略する。
【００６２】
第１特徴ベクトル算出部２３は、第１画像について特徴ベクトルを算出すると、算出した特徴ベクトルに基づいて多次元インデックスを作成する（ステップＳ７０４）。
【００６３】
式（４８）、（４９）に示すように、式（３７）、（３８）で求められる特徴ベクトルξx、ξyのノルムは１になる。
【数４８】
【数４９】
つまり、式（５０）に示すように、特徴ベクトルξxとξyの内積は、特徴ベクトルξxとξyのユークリッド距離に変換することができる。
【数５０】
つまり、特徴ベクトルξxとξyのユークリッド距離が式（５１）の条件を満たす場合、特徴ベクトルξxとξyの内積は閾値θ1以上となり、ＮＣＣの上限値ＵＢは閾値θ1以上となる。
【数５１】
【００６４】
従って、画像照合装置１は、特徴ベクトルξxとξyの内積を算出する代わりに、特徴ベクトルξxとξyのユークリッド空間における距離探索を行うことにより、式（５２）の条件を満たす画像のペアを抽出することができる。
【数５２】
【００６５】
この場合、画像照合装置１は、ユークリッド距離を用いた距離探索をサポートする、多次元データに対するインデキシング技術（Bohm, C., Berchtold, S., Keim, D. A., September 2001. Searching in high-dimensional spaces: Index structures for improving the performance of multimedia databases. ACM Computing Surveys 33 (3), 322-373.）を適用することができ、画像照合のさらなる高速化を図ることができる。画像照合装置１は、多次元インデキシング技術として、例えば、木構造を用いるＡＮＮ（Approximate Nearest Neighbors）（Arya, S., Mount, D. M., Netanyahu, N. S., Silverman, R., Wu, A. Y., November 1998. An optimal algorithm for approximate nearest neighbor searching fixed dimensions. J. ACM 45, 891-923.）、ハッシュを用いるＬＳＨ（Locality Sensitive Hashing）（Andoni, A., Indyk, P., 2008. Near-optimal hashing algorithms for approximate nearest neighbor in high dimensions. Commun. ACM 51 (1), 117-122.）、階層的ベクトル量子化（Nister, D., Stewenius, H., 2006. Scalable recognition with a vocabulary tree. In: Proc. of CVPR. Vol. 2. pp. 2161-2168.）、ベクトル量子化及びスカラー量子化（ハミングエンベディング）（Douze, M., Je´gou, H., Sandhawalia, H., Amsaleg, L., Schmid, C., 2009. Evaluation of gist descriptors for web-scale image search. In: Proceeding of the ACM International Conference on Image and Video Retrieval. CIVR ’09. ACM, New York, NY, USA, pp. 19:1-19:8.）等の、高次元空間におけるｋ近傍法（ｋ−ＮＮ探索）を用いた技術を利用できる。
【００６６】
例えば、多次元インデキシング技術としてＬＳＨを利用する場合、第１特徴ベクトル算出部２３は、特徴ベクトルに対して、予め定められた複数のハッシュ関数を用いて複数のハッシュ値を算出し、多次元インデックスとする。
【００６７】
あるいは、多次元インデキシング技術としてＡＮＮを利用する場合、画像照合装置１は、特徴ベクトルについての特徴空間を予め所定の領域に分割しておく。そして、第１特徴ベクトル算出部２３は、特徴ベクトルが属する特徴空間の領域を多次元インデックスとする。
【００６８】
第１特徴ベクトル算出部２３は、多次元インデックスを作成すると、第１画像と関連付けて記憶部１２に格納し（ステップＳ７０５）、一連のステップを終了する。
【００６９】
図８は、画像照合装置１が画像検索サーバとして動作する場合の画像照合処理の動作を示すフローチャートである。以下、図８に示したフローチャートを参照しつつ、画像照合処理の動作を説明する。なお、以下に説明する動作のフローは、予め記憶部１２に記憶されているプログラムに基づき主に制御部４０により画像照合装置１の各要素と協同して実行される。
【００７０】
なお、図８に示すステップＳ８０１〜Ｓ８０３の処理は、図２に示すステップＳ２０４〜Ｓ２０６の処理と同じであるため、説明を省略する。ただし、多次元インデキシング技術としてＬＳＨを利用する場合、第２特徴ベクトル算出部２６は、ステップＳ８０３で第２画像について特徴ベクトルを算出した後、その特徴ベクトルについて、第１画像の特徴ベクトルに対してハッシュ値を算出したのと同一の、複数のハッシュ関数を用いて複数のハッシュ値を算出する。
【００７１】
ステップＳ８０４〜Ｓ８０８の処理は、記憶部１２に格納された全ての第１画像に対して各画像毎に実施される。
【００７２】
ステップＳ８０４において、インデックス判定部４１は、第２特徴ベクトル算出部２６により第２画像について特徴ベクトルが算出されると、記憶部１２から第１画像の特徴ベクトル及び多次元インデックスを読み出す（ステップＳ８０４）。
【００７３】
次に、インデックス判定部４１は、第１画像の特徴ベクトルと第２画像の特徴ベクトルが多次元インデックスに基づく条件を満たしているか否かを判定する（ステップＳ８０５）。
【００７４】
例えば、多次元インデキシング技術としてＬＳＨを利用する場合、インデックス判定部４１は、記憶部１２から読み出した第１画像の特徴ベクトルのハッシュ値が、第２画像の特徴ベクトルの対応するハッシュ値と所定数（例えば、全ハッシュ値のうちの半数）以上一致するか否かにより、多次元インデックスに基づく条件を満たしているか否かを判定する。
【００７５】
あるいは、多次元インデキシング技術としてＡＮＮを利用する場合、インデックス判定部４１は、第１画像の特徴ベクトルが属する特徴空間の領域が、第２画像の特徴ベクトルとのユークリッド距離が（２−２θ1）以下となる部分を含むか否かにより、多次元インデックスに基づく条件を満たしているか否かを判定する。
【００７６】
インデックス判定部４１が多次元インデックスに基づく条件を満たしていると判定した場合、第１照合部２８は、第１画像の特徴ベクトルと第２画像の特徴ベクトルのユークリッド距離が式（５１）の条件を満たすか否かを判定する（ステップＳ８０６）。
【００７７】
第１照合部２８は、式（５１）の条件が満たされている場合、その第１画像を第２画像に類似する画像の候補（以下、類似画像候補と称する）として抽出する（ステップＳ８０７）。
【００７８】
一方、ステップＳ８０５でインデックス判定部４１が多次元インデックスに基づく条件が満たされていないと判定した場合、ステップＳ８０６で第１照合部２８が式（５１）の条件が満たされていないと判定した場合、又は、ステップＳ８０７で第１照合部２８が類似画像候補を抽出した場合、制御部４０は、第２画像を全ての第１画像と比較したか否かを判定する（ステップＳ８０８）。
【００７９】
制御部４０は、まだ第２画像と比較していない第１画像がある場合、ステップＳ８０４〜Ｓ８０７の処理を繰り返し、全ての第２画像と第１画像との比較が完了すると、一連のステップを終了する。
【００８０】
制御部４０は、図８のフローチャートを終了した後、さらに、図５に示したフローチャートを実施し、第２画像の類似画像候補として抽出した第１画像についてＮＣＣ値を用いて詳細に画像照合を行う。そして、制御部４０は、第２画像と類似すると判定した第１画像の情報をインターフェース部１１を介して検索を要求したコンピュータに通知する。あるいは、制御部４０は、図８のフローチャートにより類似画像候補として抽出した第１画像の情報をインターフェース部１１を介して検索を要求したコンピュータに通知してもよい。
【００８１】
以上詳述したように、図８に示したフローチャートに従って動作することによって、画像照合装置１は、画像検索サーバとして動作し、複数の画像に対する画像照合を高速かつ高精度に実施することが可能となった。また、所定の画像を複数の画像と照合する場合に多次元インデックス構造を利用して照合処理を高速化することが可能となった。
【００８２】
以下、本実施形態の低次元変換方法による類似画像の抽出精度について説明する。複数の画像Ｉi1を有する画像セットＩＳ1（ＩＳ1=｛Ｉi1｝）と複数の画像Ｉi2を有する画像セットＩＳ2（ＩＳ2=｛Ｉi2｝）において、ＮＣＣ値が閾値θ0より大きくなる画像のペアの集合ＩＩＰＳは、式（５３）のように表される。
【数５３】
この集合ＩＩＰＳは、本実施形態で算出されるＮＣＣの上限値ＵＢが閾値θ0より大きくなる、つまり式（５４）で表される画像のペアの集合ＩＩＰＳ'に近似できる。
【数５４】
【００８３】
上述した通り、ＮＣＣ値は上限値ＵＢ以下となるので、上限値ＵＢがθ0以下である場合、ＮＣＣ値は確実にθ0以下となり、対応する画像のペアは確実にＩＩＰＳに含まれないと判断できる。つまり、式（５５）の関係が成立し、ＩＩＰＳに含まれる画像のペアがＩＩＰＳ'から取りこぼされることはない。
【数５５】
【００８４】
一方、ＩＩＰＳに対してＩＩＰＳ'が大きくなりすぎると、類似しない画像まで類似画像として抽出していることになり、類似画像の抽出精度が低いということになる。このように、低次元変換の効果及び抽出した画像の類似性は、類似画像の抽出精度（以下、プレシジョンと称する）と、類似画像の非取りこぼし率（以下、リコールと称する）を用いて評価することができる。式（５６）はプレシジョンの算出式であり、式（５７）はリコールの算出式である。
【数５６】
【数５７】
例えば、ＩＩＰＳ'を非常に大きくする（極端な例ではＩＩＰＳ'をＩＩＰＳの全集合とする）ことにより、リコールを容易に１００％にすることができるが、その場合、プレシジョンが低くなる。低次元変換においてはリコールを１００％に保ちつつ、プレシジョンを可能な限り高くすることが好ましい。
【００８５】
図９は、本実施形態及び他の低次元変換方法においてリコールを１００％にしたときのプレシジョンを表すグラフである。図９では、他の低次元変換方法としてＭＳＥＡと、離散コサイン変換（ＤＣＴ：Discrete Cosine Transform）が用いられる。
【００８６】
図９に示すグラフ９００において、グラフ９０１は本実施形態によるプレシジョンを示し、グラフ９０２はＭＳＥＡによるプレシジョンを示し、グラフ９０３はＤＣＴによるプレシジョンを示す。グラフ９００の横軸は各低次元変換における特徴量（特徴ベクトル）の次元数を示し、縦軸はプレシジョンを示す。なお、図９は、１時間のビデオからランダムに選択された１万フレームの画像と他の１時間のビデオからランダムに選択された１万フレームの画像についてのプレシジョンを示す。この二つの異なる１時間のビデオは、異なった日の同じ時間帯に同じチャンネルで放送された映像であり、オープニング、エンディングの画像等、同一の画像が幾つか含まれている。各画像は３５２×２４０画素の画像であり、閾値θ0及びθ1は０．９に設定されている。
【００８７】
式（５）の上限値ＵＢMSEAは、式（５８）に示すように、式（５９）、（６０）で示される特徴ベクトルの内積で表される。
【数５８】
【数５９】
【数６０】
そこで、ＭＳＥＡを用いた低次元変換方法では、低次元特徴量として式（５９）、（６０）に示す特徴ベクトルを使用し、その内積が閾値θ1以上となる二画像を類似画像として抽出する。
【００８８】
図１０は、ＤＣＴを用いた低次元変換方法を説明するための概略図である。ＤＣＴを用いた低次元変換方法では、まず、各画像の正規化画像を所定のブロックに分割し、二次元ＤＣＴにより周波数成分に変換する。図１０のブロック１０００内の各要素（ａ0、ａ1、ａ2、...）は、二次元ＤＣＴにより変換されたＤＣＴ係数を示す。図１０に示すように、二次元ＤＣＴにより変換されたＤＣＴ係数は、ＭＰＥＧ、ＪＰＥＧ等で用いられるジグザグスキャンにより低周波数成分側から順に並べられる。そして、式（６１）に示すように、低周波数成分側から順に選択されたＤＣＴ係数による特徴ベクトルがＤＣＴによる低次元特徴量となる。
【数６１】
なお、ＤＣＴ係数は正規化画像から生成され、ＤＣ成分は常に０となるため省略される。
【００８９】
ＤＣＴにより変換された周波数成分の全パワーは元の画像のパワーと一致する。つまり、全ＤＣＴ係数を用いて計算したユークリッド距離がＮＣＣ値に対応するので、ＤＣＴ係数の一部（又は全部）を用いて計算したユークリッド距離が（２−２θ0）以下である場合、ＮＣＣ値がθ0以上となることはない。そこで、この方法では、式（６２）に示すように、それぞれのＤＣＴ成分がａi、ｂiである二つの画像の特徴ベクトルのユークリッド距離を求め、ユークリッド距離が（２−２θ1）（θ1≧θ0）以下となる二画像を類似画像として抽出する。
【数６２】
【００９０】
なお、図９に示すグラフ９００では、式（６２）に示すユークリッド距離が（２−２θ1）以下となる画像のペアの集合をＩＩＰＳ'とし、ＮＣＣ値がθ0以上となる画像のペアの集合をＩＩＰＳとしてプレシジョン及びリコールを算出している。
【００９１】
図９に示すように、（特に次元数が低い場合）プレシジョンは必ずしも高くない。プレシジョンは、最も高くても１５％未満であり（次元数が１０２４の場合）、次元数が６４の場合は５〜６％まで落ちる。また、ＤＣＴによるプレシジョンは、本実施形態によるプレシジョンとほとんど同じ値であるが、ＭＳＥＡによるプレシジョンは、極めて低い値となる。
【００９２】
一方、例えば、ＮＣＣの上限値ＵＢについての閾値θ1をＮＣＣ値についての閾値θ0より大きい値にすると（つまり、類似画像を正確にフィルタリングするのではなく、近似的にフィルタリングすると）、リコールは１００％にならなくなるが、プレシジョンを高くすることができる。
【００９３】
図１１（ａ）は、本実施形態及び他の低次元変換方法においてリコールを１００％以下にしたときのリコールとプレシジョンの関係のグラフを示し、図１１（ｂ）は、その拡大図を示す。図１１（ａ）、（ｂ）に示すグラフ１１００、１１１０では、他の低次元変換方法として、ＭＳＥＡと、ＤＣＴと、正規化輝度ヒストグラム（以下、ＮＩＨ（Normalized intensity histograms）と称する）が用いられる。
【００９４】
図１２は、ＮＩＨを用いた低次元変換方法を説明するための概略図である。ＮＩＨを用いた低次元変換方法では、まず、各画像の正規化画像１２００を（２×２、３×３等の）サブ領域１２１０に分割し、各サブ領域について、一定範囲毎に量子化した正規化輝度のヒストグラム１２２０を算出する。そして、式（６３）に示すように、全てのサブ領域のヒストグラムの分布値ａi,j（ｉはサブ領域の番号、ｊは各ヒストグラムの量子化レベルの番号）を順番に並べたベクトルがＮＩＨによる低次元特徴量となる。
【数６３】
この方法では、式（６４）に示すように、各ヒストグラムの分布値がａi,j、ｂi,jである二つの画像の特徴ベクトルのユークリッド距離を求め、ユークリッド距離が（２−２θ1）以下となる二画像を類似画像として抽出する。
【数６４】
【００９５】
また、この方法では、サブ領域の数と正規化輝度の量子化レベルを変更することにより、特徴量の次元数を変更することができる。例えば、サブ領域を２×２にし、正規化輝度の量子化を４レベルにした場合、特徴量の次元数は、２×２×４となる。ＮＩＨを用いて類似画像を抽出する場合、ＮＣＣの上限値は算出されないのでＮＣＣ値が閾値θ0以上になる画像のペアの集合ＩＩＰＳに対するリコールは１００％とならないが、リコール及びプレシジョンが非常に高くなることが知られている。
【００９６】
なお、図１１に示すグラフ１１００では、式（６４）に示すユークリッド距離が（２−２θ1）以下となる画像のペアの集合をＩＩＰＳ'とし、ＮＣＣ値がθ0以上の画像のペアの集合をＩＩＰＳとしてプレシジョン及びリコールを算出している。
【００９７】
図１１（ａ）、（ｂ）に示すグラフ１１００において、「ｐｒｏｐ １６Ｄ」、「ｐｒｏｐ ６４Ｄ」、「ｐｒｏｐ ２５６Ｄ」、「ｐｒｏｐ １０２４Ｄ」のグラフは、本実施形態の低次元変換方法で特徴量次元数がそれぞれ１６、６４、２５６、１０２４のときのリコールとプレシジョンの関係を示す。また、「ＤＣＴ ５Ｄ」、「ＤＣＴ ４４Ｄ」、「ＤＣＴ ２０９Ｄ」のグラフは、ＤＣＴによる低次元変換方法で特徴量次元数がそれぞれ５、４４、２０９のときのリコールとプレシジョンの関係を示す。また、「ＮＩＨ ２×２×４」、「ＮＩＨ ４×４×４」、「ＮＩＨ ４×４×８」のグラフは、ＮＩＨによる低次元変換方法で特徴量次元数がそれぞれ２×２×４、４×４×４、４×４×８のときのリコールとプレシジョンの関係を示す。また、「ＭＳＥＡ １６Ｄ」、「ＭＳＥＡ ６４Ｄ」、「ＭＳＥＡ ５１２Ｄ」のグラフは、ＭＳＥＡによる低次元変換方法で特徴量次元数がそれぞれ１６、６４、５１２のときのリコールとプレシジョンの関係を示す。グラフ１１００、１１１０の横軸はリコールの値を示し、縦軸はプレシジョンの値を示す。
【００９８】
図１１（ａ）、（ｂ）に示すように、本実施形態の低次元変換方法では、特徴量次元数が６４と２５６の場合、リコールを９７％に保ちつつ、プレシジョンが７０％以上となる。特に、リコールが非常に高い（９５％以上）場合、ＤＣＴ、ＮＩＨ、ＭＳＥＡよりプレシジョンが高くなる。
【００９９】
図１３（ａ）は、図１１（ａ）とは異なる画像セットのペア（それぞれ、類似する放送映像から取得した１万画像）についてのリコールとプレシジョンの関係のグラフを示し、図１３（ｂ）は、その拡大図を示す。図１３（ａ）、（ｂ）のグラフ１３００、１３１０では、ほとんど全てのタイプで、図１１（ａ）、（ｂ）のグラフ１１００、１１１０より高いリコールとプレシジョンの組合せが得られる。図１３（ｂ）に示すように、本実施形態の低次元変換方法で特徴量次元数が１０２４又は２５６のときに最も高いパフォーマンスが得られ、次いでＭＳＥＡによる低次元変換方法で特徴量次元数が５１２のとき、本実施形態の低次元変換方法で特徴量次元数が６４又は１６のとき、ＭＳＥＡによる低次元変換方法で特徴量次元数が６４のときの順に高いパフォーマンスが得られる。一方、ＤＣＴによる低次元変換方法では、あまり高いパフォーマンスが得られない。また、ＭＳＥＡによる低次元変換方法で特徴量次元数が１６の場合のパフォーマンスは、本実施形態の低次元変換方法で特徴量次元数が１６の場合より非常に低くなる。つまり、本実施形態の低次元変換方法は、全ての特徴量次元数において、他の低次元変換方法より高いパフォーマンスを得ることができる。
【０１００】
図１４（ａ）、（ｂ）、（ｃ）、（ｄ）は、それぞれリコールを０．９５、０．９７、０．９８、０．９９に固定したときの各低次元変換方法での特徴量次元数とプレシジョンの関係のグラフを示す。図１４（ａ）、（ｂ）、（ｃ）、（ｄ）に示すグラフ１４００、１４１０、１４２０、１４３０において、それぞれ、グラフ１４０１、１４１１、１４２１、１４３１は本実施形態の低次元変換方法でのプレシジョンを示し、グラフ１４０２、１４１２、１４２２、１４３２はＭＳＥＡによる低次元変換方法でのプレシジョンを示し、グラフ１４０３、１４１３、１４２３、１４３３はＤＣＴによる低次元変換方法でのプレシジョンを示し、グラフ１４０４、１４１４、１４２４、１４３４はＮＩＨによる低次元変換方法でのプレシジョンを示す。グラフ１４００〜１４３０の横軸は特徴量次元数を示し、縦軸はプレシジョンを示す。
【０１０１】
図１４（ａ）〜（ｄ）に示すように、リコールを非常に高い値にした場合、本実施形態の低次元変換方法でのプレシジョンは、他の低次元変換方法でのプレシジョンより高くなることが明らかである。図１４（ｄ）に示すように、リコールが０．９９の場合、本実施形態の低次元変換方法では、ＤＣＴによる低次元変換方法よりわずかにパフォーマンスが低くなるが、この場合、そもそもプレシジョンが０．１４以下であり、極めて低い。また、例えば、図１４（ａ）に示すように、リコールが０．９５の場合、本実施形態の低次元変換方法において、特徴量次元数が低い方が高いパフォーマンスを示している（特徴量次元数が６４、１２８、２５６の方が、特徴量次元数が５１２のときよりもプレシジョンが高い）。類似した現象がＭＳＥＡでもみられる（特徴量次元数が１６、３２、６４の方が、特徴量次元数が１２８のときよりもプレシジョンが高い）。これは、これらの方法が複数ブロックへ分割するものであるためと考えられる。そのため、リコールを１００％としない場合、適切なサイズのブロックに分割する必要がある。図１４（ａ）〜（ｄ）から、本実施形態の低次元変換方法でのプレシジョンは、特徴量次元数が６４のとき、すなわち、分割するブロック数が３２のとき、最も高くなると考えられる。
【０１０２】
図１５（ａ）〜（ｄ）は、図１４（ａ）〜（ｄ）とは異なる画像セットのペアについて、それぞれリコールを０．９５、０．９７、０．９８、０．９９に固定したときの各低次元変換方法での特徴量次元数とプレシジョンの関係のグラフを示す。図１５（ａ）、（ｂ）、（ｃ）、（ｄ）に示すグラフ１５００、１５１０、１５２０、１５３０において、それぞれ、グラフ１５０１、１５１１、１５２１、１５３１は本実施形態の低次元変換方法でのプレシジョンを示し、グラフ１５０２、１５１２、１５２２、１５３２はＭＳＥＡによる低次元変換方法でのプレシジョンを示し、グラフ１５０３、１５１３、１５２３、１５３３はＤＣＴによる低次元変換方法でのプレシジョンを示し、グラフ１５０４、１５１４、１５２４、１５３４はＮＩＨによる低次元変換方法でのプレシジョンを示す。グラフ１５００〜１５３０の横軸は特徴量次元数を示し、縦軸はプレシジョンを示す。
【０１０３】
図１５（ａ）〜（ｄ）に示すグラフ１５００〜１５３０では、プレシジョンが比較的高くなっている。例えば、図１５（ｄ）では、リコールが０．９９に固定されているにも関わらず、プレシジョンは略０．６となる。図１５（ａ）〜（ｃ）に示すように、リコールが比較的低い（０．９５〜０．９８）場合、本実施形態の低次元変換方法とＭＳＥＡによる低次元変換方法はＤＣＴ及びＮＩＨによる低次元変換方法よりパフォーマンスがよい。一方、図１５（ｄ）に示すように、リコールが高い（０．９９）場合、本実施形態の低次元変換方法とＤＣＴ及びＮＩＨによる低次元変換方法はＭＳＥＡによる低次元変換方法よりパフォーマンスが高い。つまり、全体として、本実施形態の低次元変換方法は他の低次元変換方法より高いパフォーマンスを得ることができる。
【０１０４】
以上、図９、図１１（ａ）、（ｂ）、図１３（ａ）、（ｂ）、図１４（ａ）〜（ｄ）、図１５（ａ）〜（ｄ）を用いて説明したように、本実施形態の低次元変換方法を用いることにより、他の低次元変換方法より高精度に類似性の高い画像を抽出することが可能となる。
【０１０５】
本実施形態の低次元変換方法では、ラグランジュの未定乗数決定法を用いてＮＣＣの上限値を算出する。これにより、ＮＣＣ値が上限値に近いとき、そのＮＣＣの微分係数は非常に０に近くなり、上限値から離れているとき、微分係数は０から離れる。そのため、ＮＣＣ値が上限値に近い（リコールが高い）場合、上限値をわずかに大きくしても、微分係数が０に近いのでリコールは大きく変化しない。一方、ＮＣＣ値が上限値から離れている場合、微分係数は０から離れているので、プレシジョンは大幅に増加する。これにより、本実施形態の低次元変換方法では、リコールを高く保ちつつ、プレシジョンを高くできると考えられる。
【０１０６】
図１６は、本実施形態による照合処理の時間を示す表である。図１６の表１６００は、それぞれ１万画像（３５２×２４０画素）からなる二つの画像セットにおいて、全ての画像ペアの組合せについて特徴ベクトルのユークリッド距離により全数照合したときの時間を示す。なお、この時間には、特徴ベクトルの算出時間は含まれない。
【０１０７】
表１６００の行１６１０と行１６３０は低次元特徴量による照合処理を、行１６２０と行１６４０は低次元特徴量による照合処理及びそれにより抽出した類似画像のペアに対するＮＣＣによる照合処理を示す。また、行１６１０と行１６２０ではリコールが１００％になるようにし（低次元特徴量による照合処理における閾値θ1＝ＮＣＣによる照合処理における閾値θ0＝０．９）、行１６３０と行１６４０ではリコールが９５％になるようにしている（θ1＞θ0）。各行において、さらに特徴量の次元数毎に処理結果が示される。一方、列１６５０は特徴量次元数を、列１６６０はプレシジョンを、列１６７０はリコールを、列１６８０は処理時間（秒）をそれぞれ示す。
【０１０８】
行１６１０と行１６３０に示すように、照合処理時間は、概ね特徴量の次元数に比例している。これに基づくと、全ての画像ペアの組合せについてＮＣＣによる照合処理（つまり、３５２×２４０次元の演算）を実施した場合、約４０時間かかると考えられる。リコールが１００％の場合、行１６１０に示すようにプレシジョン値は比較的低くなり、類似画像として多くの画像が抽出されるため、行１６２０に示すようにＮＣＣによる照合処理時間は非常に長くなる。一方、リコールが９５％の場合、行１６３０に示すようにプレシジョンが大幅に高くなるため、類似画像の抽出数が抑制され、行１６４０に示すようにＮＣＣによる照合処理時間はあまり長くならない。従って、わずかにリコールを減少させる（９５％）のがリーズナブルであると考えられる。
【０１０９】
なお、上述したように、行１６１０と行１６３０の照合処理時間は、特徴量の次元数に応じて異なるものであり、本実施形態とＭＳＥＡとで差は生じない。また、１万画像（３５２×２４０画素）について、本実施形態に基づく低次元特徴量の算出時間とＭＳＥＡに基づく低次元特徴量の算出時間は、その次元数に関わらず、それぞれ、約３５秒と約３２秒であり、ほとんど差がなかった。つまり、次元数が同じ場合、本実施形態の低次元変換方法とＭＳＥＡによる低次元変換方法とで類似画像を抽出するのにかかる時間は、ほとんど同じとなる。一方、本実施形態の低次元変換方法では、ＭＳＥＡによる低次元変換方法より、類似画像として抽出する画像のペアが少なくなるので、その後のＮＣＣによる照合処理時間は短くなり、トータルの照合時間は短くなる。
【０１１０】
図１７は、画像照合装置を画像符号化装置に適用する場合の例を示す概略構成図である。図１７に示す画像照合装置２は、図１に示す画像照合装置１の各部に加えて画像入力部５３を有する。また、画像照合装置２の制御部６０は、インテグラル画像生成部６１、符号化処理部６２、画像分割部６３、特徴ベクトル算出部６４、上限値算出部６５、第１照合部６６、ＮＣＣ値算出部６７及び第２照合部６８を有する。なお、画像分割部６３は、図１に示す第１画像分割部２２と第２画像分割部２５の機能を備え、特徴ベクトル算出部６４は、図１に示す第１特徴ベクトル算出部２３と第２特徴ベクトル算出部２６の機能を備える。
【０１１１】
画像入力部５３は、ＣＣＤ、ＣＭＯＳ等の光電変換器で構成された２次元検出器と、その２次元検出器上に撮影対象の像を結像する結像光学系等を有する。画像入力部５３は、一定の時間間隔（例えば１／３０秒）毎に撮影を行い、撮影画像を、例えば３５２×２４０画素のデジタル画像に変換し、そのデジタル画像を記憶部５２に記憶する。
【０１１２】
図１８は、図１７に示す画像照合装置２による画像の符号化処理の動作を示すフローチャートである。以下、図１８に示したフローチャートを参照しつつ、画像符号化処理の動作を説明する。なお、以下に説明する動作のフローは、予め記憶部５２に記憶されているプログラムに基づき主に制御部６０により画像照合装置２の各要素と協同して実行される。
【０１１３】
最初に、画像入力部５３は、撮影対象を撮影した画像をデジタル画像に変換し、そのデジタル画像を記憶部５２に記憶する（ステップＳ１８０１）。
【０１１４】
次に、インテグラル画像生成部６１は、記憶部１２に保存されたデジタル画像を読み出し、インテグラル画像を生成し、記憶部５２に記憶する（ステップＳ１８０２）。
【０１１５】
デジタル画像及びインテグラル画像の水平方向をｘ軸、垂直方向をｙ軸とし、デジタル画像の座標（ａ、ｂ）における正規化画素値をＩ（ａ、ｂ）とする。インテグラル画像生成部６７は、各画素が式（６５）、（６６）からなるインテグラル画像（Viola, P., Jones, M., May 2004. Robust real-time face detection. International Journal of Computer Vision 57 (2), 137-154.）をそれぞれ生成する。
【数６５】
【数６６】
つまり、式（６５）により算出されるインテグラル画像の各画素値は、原点からその画素までの全ての画素の正規化画素値の総和であり、式（６６）により算出されるインテグラル画像の各画素値は、原点からその画素までの全ての画素の正規化画素値の二乗和である。
【０１１６】
図１９は、インテグラル画像を説明するための概略図である。例えば、インテグラル画像１９００において、座標（ｘ1、ｙ1）、（ｘ1、ｙ2）、（ｘ2、ｙ1）、（ｘ2、ｙ2）で囲まれた領域Ｄの正規化画素値の総和と正規化画素値の二乗和は、それぞれ、式（６７）、（６８）により算出することができる。
【数６７】
【数６８】
従って、予めインテグラル画像を作成しておくことにより、所定領域内の正規化画素値の総和と二乗和の算出処理を高速化できるので、特徴ベクトルの算出を高速化できる。
【０１１７】
次に、符号化処理部６２は、記憶部１２に保存されたデジタル画像を読み出し、デジタル画像のフォーマット（画素数）変換、８×８画素、１６×１６画素等の符号化ブロックへの分割等の符号化前処理を実施する（ステップＳ１８０３）。
【０１１８】
以下のステップＳ１８０４〜Ｓ１８１４の処理は、符号化ブロック毎に実施される。まず、画像分割部６３は、符号化ブロックをさらに所定サイズの複数のブロックに分割する（ステップＳ１８０４）。例えば、画像分割部６３は、１６×１６画素の画像を４×４画素の１６ブロックに分割する。
【０１１９】
次に、特徴ベクトル算出部６４は、インテグラル画像を用いて、符号化ブロックについて、式（３７）に示す特徴ベクトルを算出し、符号化ブロックと関連付けて記憶部１２に保存する（ステップＳ１８０５）。
【０１２０】
次に、上限値算出部６５は、記憶部１２から照合元ブロックの特徴ベクトルを取得する（ステップＳ１８０６）。
【０１２１】
この照合元ブロックは、例えば、１フレーム前のデジタル画像における各符号化ブロックである。なお、前方向予測だけでなく、後方向予測も用いる場合、上限値算出部６５は、１フレーム後のデジタル画像における各符号化ブロックも照合元ブロックとして用いる。
【０１２２】
次に、上限値算出部６５は、符号化ブロックの特徴ベクトルと照合元ブロックの特徴ベクトルの内積を算出して、符号化ブロックと照合元ブロックについて、ＮＣＣの上限値ＵＢを求める（ステップＳ１８０７）。
【０１２３】
次に、第１照合部６６は、符号化ブロックと照合元ブロックについての上限値ＵＢが閾値θ1以上であるか否かを判定する（ステップＳ１８０８）。
【０１２４】
上限値ＵＢが閾値θ1以上である場合、第１照合部６５は、その照合元ブロックを符号化ブロックと類似すると判定し、ＮＣＣ値算出部６７は、その符号化ブロックと照合元ブロックについてＮＣＣ値を算出する（ステップＳ１８０９）。
【０１２５】
次に、第２照合部６８は、ＮＣＣ値算出部６７が算出したＮＣＣ値がその符号化ブロックについて算出されたＮＣＣ値のうち最大であるか否かを判定する（ステップＳ１８１０）。
【０１２６】
ＮＣＣ値算出部６７が算出したＮＣＣ値がその符号化ブロックについて算出されたＮＣＣ値のうち最大である場合、第２照合部６８は、その照合元ブロックを動き検出用ブロックとして記憶部５２に記憶（すでに記憶されている場合、更新）する（ステップＳ１８１１）。
【０１２７】
一方、上限値ＵＢが閾値θ1未満である場合、ＮＣＣ値算出部６７が算出したＮＣＣ値がその符号化ブロックについて算出されたＮＣＣ値のうち最大でない場合、又は、第２照合部６８が照合元ブロックを動き検出用ブロックとして記憶部５２に記憶した場合、制御部６０は、符号化ブロックを全ての照合元ブロックと比較したか否かを判定する（ステップＳ１８１２）。
【０１２８】
制御部６０は、まだ比較していない照合元ブロックがある場合、ステップＳ１８０７〜Ｓ１８１１の処理を繰り返し、全ての照合元ブロックとの比較が完了すると、符号化処理部６２は、記憶部５２に記憶された動き検出用ブロックの候補を用いて動き補償を行う。また、符号化処理部６２は、ＤＣＴ変換、量子化等の符号化処理を実施する（ステップＳ１８１３）。
【０１２９】
次に、制御部６０は、全ての符号化ブロックの符号化処理が完了したか否かを判定する（ステップＳ１８１４）。全ての符号化ブロックの符号化処理が完了していない場合、ステップＳ１８０３〜Ｓ１８１３の処理を繰り返し、全ての符号化ブロックの符号化処理が完了すると、制御部６０は、一連のステップを終了する。
【０１３０】
このようにして符号化されたデータは、インターフェース部５１を介して、外部のデコード装置（不図示）に送信される。
【０１３１】
なお、第２照合部６８は、符号化ブロックと全ての照合元ブロックについての上限値ＵＢが閾値θ1未満であった場合は、その符号化ブロックに類似する照合元ブロックが存在しないと判断し、任意の照合元ブロックを動き検出用ブロックとして選択する。あるいは、その場合、上限値ＵＢが最も高かった照合元ブロックを動き検出用ブロックとして選択してもよい。
【０１３２】
以上詳述したように、図１８に示したフローチャートに従って動作することによって、画像照合装置２は、画像符号化装置として動作し、高速に動き検出処理を実施することが可能となった。また、インテグラル画像を用いて特徴量算出処理を高速化することが可能となった。また、本実施形態による方法は、スライディングウィンドウに基づく探索（Wei, S.-D., Lai, S.-H., 2007. Efficient normalized cross correlation based on adaptive multilevel successive elimination. In: Proceedings of the 8th Asian conference on Computer vision - Volume Part I. ACCV’07. Springer-Verlag, Berlin, Heidelberg, pp. 638-646.）とも親和性が高い。
【符号の説明】
【０１３３】
１、２ 画像照合装置
１１、５１ インターフェース部
１２、５２ 記憶部
２０、３０、４０、６０ 制御部
２１ 第１画像取得部
２２ 第１画像分割部
２３ 第１特徴ベクトル算出部
２４ 第２画像取得部
２５ 第２画像分割部
２６ 第２特徴ベクトル算出部
２７、６５ 上限値算出部
２８、６６ 第１照合部
３１、６７ ＮＣＣ値算出部
３２、６８ 第２照合部
４１ インデックス判定部
５３ 画像入力部
６１ インテグラル画像生成部
６２ 符号化処理部
６３ 画像分割部
６４ 特徴ベクトル算出部
【特許請求の範囲】
【請求項１】
第１画像を複数のブロックに分割し、第２画像を前記第１画像と同数のブロックに分割する画像分割部と、
前記第１画像及び前記第２画像のそれぞれについて、各ブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づき、前記第１画像の特徴ベクトル及び前記第２画像の特徴ベクトルを算出する特徴ベクトル算出部と、
ラグランジュの未定乗算決定法を利用して２つの画像の正規化相互相関値を求める式から導いた２つの画像の正規化相互相関の上限値を算出する式を用いて、算出された前記第１画像の特徴ベクトル及び前記第２画像の特徴ベクトルに基づき、前記第１画像と前記第２画像の正規化相互相関の上限値を算出する上限値算出部と、
前記上限値が第１の閾値以上であるか否かに基づいて前記第１画像と前記第２画像を照合する第１照合部と、
を有することを特徴とする画像照合装置。
【請求項２】
前記特徴ベクトル算出部は、ｍを前記ブロック内の画素数、ｎを前記ブロックの数、φi,jを前記第１画像におけるｉ番目のブロック内のｊ番目の画素の画素値、φaveを前記第１画像の全画素の画素値の平均値として、次のベクトルξxを前記第１画像の特徴ベクトルとして算出し、
【数１】
ψi,jを前記第２画像におけるｉ番目のブロック内のｊ番目の画素の画素値、ψaveを前記第２画像の全ての画素の画素値の平均値として、次のベクトルξyを前記第２画像の特徴ベクトルとして算出し、
【数２】
前記上限値算出部は、前記第１画像の特徴ベクトルと前記第２画像の特徴ベクトルの内積を前記上限値として算出する、請求項１に記載の画像照合装置。
【請求項３】
前記上限値が前記第１の閾値以上である前記第１画像と前記第２画像の正規化相互相関値を算出する正規化相互相関値算出部と、
前記正規化相互相関値が前記第１の閾値以下の値である第２の閾値以上であるか否かに基づいて前記第１画像と前記第２画像を詳細に照合する第２照合部と、をさらに有する、請求項１または２に記載の画像照合装置。
【請求項４】
第１画像を複数のブロックに分割するステップと、
前記第１画像について、各ブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づき、前記第１画像の特徴ベクトルを算出するステップと、
第２画像を前記第１画像と同数のブロックに分割するステップと、
前記第２画像について、各ブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づき、前記第２画像の特徴ベクトルを算出するステップと、
ラグランジュの未定乗算決定法を利用して２つの画像の正規化相互相関値を求める式から導いた２つの画像の正規化相互相関の上限値を算出する式を用いて、算出された前記第１画像の特徴ベクトル及び前記第２画像の特徴ベクトルに基づき、前記第１画像と前記第２画像の正規化相互相関の上限値を算出するステップと、
前記上限値が第１の閾値以上であるか否かに基づいて前記第１画像と前記第２画像を照合するステップと、
を含むことを特徴とする画像照合方法。
【請求項５】
前記第１画像の特徴ベクトルを算出するステップにおいて、ｍを前記ブロック内の画素数、ｎを前記ブロックの数、φi,jを前記第１画像におけるｉ番目のブロック内のｊ番目の画素の画素値、φaveを前記第１画像の全画素の画素値の平均値として、次のベクトルξxを前記第１画像の特徴ベクトルとして算出し、
【数３】
前記第２画像の特徴ベクトルを算出するステップにおいて、ｍを前記ブロック内の画素数、ｎを前記ブロックの数、ψi,jを前記第２画像におけるｉ番目のブロック内のｊ番目の画素の画素値、ψaveを前記第２画像の全ての画素の画素値の平均値として、次のベクトルξyを前記第２画像の特徴ベクトルとして算出し、
【数４】
前記上限値を算出するステップにおいて、前記第１画像の特徴ベクトルと前記第２画像の特徴ベクトルの内積を前記上限値として算出する、請求項４に記載の画像照合方法。
【請求項６】
前記上限値が前記第１の閾値以上である前記第１画像と前記第２画像の正規化相互相関値を算出するステップと、
前記正規化相互相関値が前記第１の閾値以下の値である第２の閾値以上であるか否かに基づいて前記第１画像と前記第２画像を詳細に照合するステップと、をさらに有する、請求項４または５に記載の画像照合方法。
【請求項７】
第１画像を複数のブロックに分割するステップと、
前記第１画像について、各ブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づき、前記第１画像の特徴ベクトルを算出するステップと、
第２画像を前記第１画像と同数のブロックに分割するステップと、
前記第２画像について、各ブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づき、前記第２画像の特徴ベクトルを算出するステップと、
ラグランジュの未定乗算決定法を利用して２つの画像の正規化相互相関値を求める式から導いた２つの画像の正規化相互相関の上限値を算出する式を用いて、算出された前記第１画像の特徴ベクトル及び前記第２画像の特徴ベクトルに基づき、前記第１画像と前記第２画像の正規化相互相関の上限値を算出するステップと、
前記上限値が第１の閾値以上であるか否かに基づいて前記第１画像と前記第２画像を照合するステップと、
をコンピュータに実行させることを特徴とするコンピュータプログラム。
【請求項１】
第１画像を複数のブロックに分割し、第２画像を前記第１画像と同数のブロックに分割する画像分割部と、
前記第１画像及び前記第２画像のそれぞれについて、各ブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づき、前記第１画像の特徴ベクトル及び前記第２画像の特徴ベクトルを算出する特徴ベクトル算出部と、
ラグランジュの未定乗算決定法を利用して２つの画像の正規化相互相関値を求める式から導いた２つの画像の正規化相互相関の上限値を算出する式を用いて、算出された前記第１画像の特徴ベクトル及び前記第２画像の特徴ベクトルに基づき、前記第１画像と前記第２画像の正規化相互相関の上限値を算出する上限値算出部と、
前記上限値が第１の閾値以上であるか否かに基づいて前記第１画像と前記第２画像を照合する第１照合部と、
を有することを特徴とする画像照合装置。
【請求項２】
前記特徴ベクトル算出部は、ｍを前記ブロック内の画素数、ｎを前記ブロックの数、φi,jを前記第１画像におけるｉ番目のブロック内のｊ番目の画素の画素値、φaveを前記第１画像の全画素の画素値の平均値として、次のベクトルξxを前記第１画像の特徴ベクトルとして算出し、
【数１】
ψi,jを前記第２画像におけるｉ番目のブロック内のｊ番目の画素の画素値、ψaveを前記第２画像の全ての画素の画素値の平均値として、次のベクトルξyを前記第２画像の特徴ベクトルとして算出し、
【数２】
前記上限値算出部は、前記第１画像の特徴ベクトルと前記第２画像の特徴ベクトルの内積を前記上限値として算出する、請求項１に記載の画像照合装置。
【請求項３】
前記上限値が前記第１の閾値以上である前記第１画像と前記第２画像の正規化相互相関値を算出する正規化相互相関値算出部と、
前記正規化相互相関値が前記第１の閾値以下の値である第２の閾値以上であるか否かに基づいて前記第１画像と前記第２画像を詳細に照合する第２照合部と、をさらに有する、請求項１または２に記載の画像照合装置。
【請求項４】
第１画像を複数のブロックに分割するステップと、
前記第１画像について、各ブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づき、前記第１画像の特徴ベクトルを算出するステップと、
第２画像を前記第１画像と同数のブロックに分割するステップと、
前記第２画像について、各ブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づき、前記第２画像の特徴ベクトルを算出するステップと、
ラグランジュの未定乗算決定法を利用して２つの画像の正規化相互相関値を求める式から導いた２つの画像の正規化相互相関の上限値を算出する式を用いて、算出された前記第１画像の特徴ベクトル及び前記第２画像の特徴ベクトルに基づき、前記第１画像と前記第２画像の正規化相互相関の上限値を算出するステップと、
前記上限値が第１の閾値以上であるか否かに基づいて前記第１画像と前記第２画像を照合するステップと、
を含むことを特徴とする画像照合方法。
【請求項５】
前記第１画像の特徴ベクトルを算出するステップにおいて、ｍを前記ブロック内の画素数、ｎを前記ブロックの数、φi,jを前記第１画像におけるｉ番目のブロック内のｊ番目の画素の画素値、φaveを前記第１画像の全画素の画素値の平均値として、次のベクトルξxを前記第１画像の特徴ベクトルとして算出し、
【数３】
前記第２画像の特徴ベクトルを算出するステップにおいて、ｍを前記ブロック内の画素数、ｎを前記ブロックの数、ψi,jを前記第２画像におけるｉ番目のブロック内のｊ番目の画素の画素値、ψaveを前記第２画像の全ての画素の画素値の平均値として、次のベクトルξyを前記第２画像の特徴ベクトルとして算出し、
【数４】
前記上限値を算出するステップにおいて、前記第１画像の特徴ベクトルと前記第２画像の特徴ベクトルの内積を前記上限値として算出する、請求項４に記載の画像照合方法。
【請求項６】
前記上限値が前記第１の閾値以上である前記第１画像と前記第２画像の正規化相互相関値を算出するステップと、
前記正規化相互相関値が前記第１の閾値以下の値である第２の閾値以上であるか否かに基づいて前記第１画像と前記第２画像を詳細に照合するステップと、をさらに有する、請求項４または５に記載の画像照合方法。
【請求項７】
第１画像を複数のブロックに分割するステップと、
前記第１画像について、各ブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づき、前記第１画像の特徴ベクトルを算出するステップと、
第２画像を前記第１画像と同数のブロックに分割するステップと、
前記第２画像について、各ブロックに含まれる画素の正規化画素値の総和及び各ブロックに含まれる画素の正規化画素値の二乗和に基づき、前記第２画像の特徴ベクトルを算出するステップと、
ラグランジュの未定乗算決定法を利用して２つの画像の正規化相互相関値を求める式から導いた２つの画像の正規化相互相関の上限値を算出する式を用いて、算出された前記第１画像の特徴ベクトル及び前記第２画像の特徴ベクトルに基づき、前記第１画像と前記第２画像の正規化相互相関の上限値を算出するステップと、
前記上限値が第１の閾値以上であるか否かに基づいて前記第１画像と前記第２画像を照合するステップと、
をコンピュータに実行させることを特徴とするコンピュータプログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【公開番号】特開２０１３−１６０７３（Ｐ２０１３−１６０７３Ａ）
【公開日】平成２５年１月２４日（２０１３．１．２４）
【国際特許分類】
【出願番号】特願２０１１−１４９４２８（Ｐ２０１１−１４９４２８）
【出願日】平成２３年７月５日（２０１１．７．５）
【出願人】（５０４２０２４７２）大学共同利用機関法人情報・システム研究機構 (119)
【Ｆターム（参考）】
【公開日】平成２５年１月２４日（２０１３．１．２４）
【国際特許分類】
【出願日】平成２３年７月５日（２０１１．７．５）
【出願人】（５０４２０２４７２）大学共同利用機関法人情報・システム研究機構 (119)
【Ｆターム（参考）】
[ Back to top ]