本願発明は、質の測定基準を計算する前に、歪画像を基準画像にアライメントする、歪画像を前処理するシステムおよび方法を提示する。基準画像を歪画像にアライメントする処理（又は、基準画像を歪画像に揃える処理、基準画像を歪画像に位置合わせする処理）は、「レジストレーションすること（registering）」または「レジストレーション（registration）」とも呼ばれることがある。

歪画像の画質評価の完全参照（ｆｕｌｌ ｒｅｆｅｒｅｎｃｅ）方法は、歪画像と基準画像との特別な形態の比較に基づく。歪みが２つの画像の間の寸法の違い、その包絡線内での歪画像の水平方向および／または垂直方向の変位または回転などの幾何学的特性の変化、または、他の変形を含む場合、既存の完全参照品質評価方法では、視覚的な評価に対応し得る品質メトリック（測定基準）を生成しない。

本発明の実施形態で提示された画像レジストレーション方法の他の応用は、視覚的観察では（ほとんど）同一に見えるが、スケール、回転、または一方の画像にノイズが多いなどの点で異なる、２つの画像の自動比較である。

したがって、画像レジストレーションはリモートセンシング（地図作製の更新）、およびコンピュータビジョンに応用される。

（異なる時点で撮影された同一の患者の画像データ、例えば、変化の検出や腫瘍の監視など、についての）医療画像レジストレーションが、例えば、呼吸または解剖学上の変化による弾性（非剛性（ｎｏｎｒｉｇｉｄ）としても知られている）変形を補償するために用いられ得る。医療画像の非剛性のレジストレーションは、解剖学上のアトラス（ａｔｌａｓ）、例えば神経画像のためのタライラッハアトラス（Ｔａｌａｉｒａｃｈ ａｔｌａｓ）などに患者の画像データをレジストレーションするためにも役に立つ。

画像レジストレーションのさらなる使用は、天体写真において宇宙を撮影した画像をアライメントすることにある。（自動または手動で入力される）制御点を用いて、一の画像上で主要な特徴を第２の画像にアライメントするように変形を行い得る。

また、画像レジストレーションはパノラマ状の画像生成の本質的な部分である。本発明に提示された技術は、カメラおよびカメラつき携帯電話などの組み込み機器においてリアルタイムで実装されて実施されている。

図１は、基準画像Ｘを提供する基準画像源１１０、歪画像Ｙを提供する歪み画像源１２０、および、アライメントパラメータ推定処理１４０および画像アライメント処理１５０を含む画像前処理システム１３０、および、画像評価処理１６０を有する画像アライメントシステム１００の概略ブロック図を示す。

基準画像Ｘおよび歪画像Ｙの双方が、画像前処理システム１３０に入力される。ここで、歪画像Ｙは基準画像Ｘを参照して前処理され、前処理された画像Ｙ^＊を生成し、それはレジストレーション画像Ｙ^＊とも呼ばれる。

概して、画像評価処理１６０は入力として、基準画像Ｘおよび前処理画像Ｙ^＊を受信し、評価結果Ｑ０を生成する。

歪画像Ｙは、基準画像Ｘにさかのぼることができる場合があり、歪画像Ｙの歪みは、圧縮と、その後の画像の圧縮解除；画像の回転、平行移動、または線形スケーリングによって特徴づけられる幾何学的変形；画像のクロッピングによる画像の一部の喪失；ノイズのあるチャネル上での転送によるノイズの追加のいずれかの単独または組み合わせの結果である場合がある。

アライメントパラメータ推定処理１４０の目的は、基準画像Ｘを参照して歪画像Ｙを解析し、歪みの幾何学的特徴の推定を提供するアライメントパラメータＡＰの組を生成することである。

アライメントパラメータＡＰは歪画像Ｙを前処理された画像Ｙ^＊に変形する画像アライメント処理１５０によって使用される。前処理された画像Ｙ^＊はレジストレーション画像Ｙ^＊とも呼ばれることがある。アライメントパラメータＡＰを用いて、画像アライメントプロセス１５０は歪画像Ｙのいくらかまたはすべての幾何学的歪みを補償するが、画像それ自体の視覚的性質は変更しない。結果として、歪画像Ｙは画像アライメント処理１５０によってレジストレーション画像Ｙ^＊に変形される。すなわち、それは、基準画像Ｘにレジストレーションされてアライメントされる。

画像評価処理１６０は、従来の視覚品質測定に基づいていてもよく、その場合、評価結果Ｑ０は、構造的類似性（ＳＳＩＭ）指標または視覚情報忠実度（ＶＩＦ）指標などの品質の測定基準（品質メトリック）を表す。この場合、評価結果Ｑ０は、前処理画像Ｙ^＊に関して計算されるが、歪画像Ｙの真の視覚的品質メトリックとして仮定される。なぜなら、前処理の行為は画像の視覚的内容を変更することなく幾何学的歪みを補正するだけであるためである。例えば、ノイズ及び圧縮の効果は残る。

多くの他の応用において、画像評価処理１６０の出力Ｑ０は、単に、歪画像Ｙと基準画像Ｘの間の、例えば両方の画像が同一の場面の別々の写真である場合の、視覚的類似性の指標であり得る。

画像前処理システム１３０の第１の実施形態は、以下に基づく。
−基準画像Ｘと歪画像Ｙの両方からキーポイント（keypoint）の組を抽出するステップ。
−最も近い、ＸおよびＹから選択された数のキーポイントのペアをマッチングさせるステップ。
−第１のアフィン変換のパラメータに関して解くステップ。
−第１のアフィン変換の逆（ｉｎｖｅｒｓｅ）である第２のアフィン変換を生成するステップ。
−歪画像Ｙに第２のアフィン変換を適用し、それによって前処理された画像Ｙ^＊を生成するステップ。

画像内の重要な領域を特徴づけるキーポイントは、Ｄ．Ｇ．Ｌｏｗｅ「Ｄｉｓｔｉｎｃｔｉｖｅ Ｉｍａｇｅ Ｆｅａｔｕｒｅｓ ｆｒｏｍ Ｓｃａｌｅ−Ｉｎｖａｒｉａｎｔ Ｋｅｙｐｏｉｎｔｓ」、Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｊｏｕｒｎａｌ ｏｆ Ｃｏｍｐｕｔｅｒ Ｖｉｓｉｏｎ、第６０巻、第２号、第９１−１１０ページ、２００４年によって提示されたスケール不変特徴変換（ＳＩＦＴ）パラメータセットとして知られている。アライメントパラメータ推定処理１４０の２つの特定の実施形態および画像アライメント処理１５０の対応する実施形態が提示される。

図２ａは、画像前処理システム１３０の第１の実施形態である、アフィン変換に基づく画像アライメントシステム２００．１のブロック図を示す。システム２００．１は、図１のアライメントパラメータ推定処理１４０の実施形態であるアフィン変換パラメータ推定サブシステム２１０および、図１の画像アライメント処理１５０の実施形態であるアフィン画像変換モジュール２２０を含む。

アフィン変換パラメータ推定サブシステム２１０は、画素埋込モジュール２３０、Ｘ−ＳＩＦＴ抽出モジュール２４０、Ｙ−ＳＩＦＴ抽出モジュール２５０、キーポイントマッチングモジュール２６０、アフィン変換生成モジュール２７０、および変換反転モジュール２８０を含む。

画素埋込モジュール２３０においては、画素で計測された水平方向寸法「Ｍ」および垂直方向寸法「Ｎ」に関して２つの画像ＸおよびＹが同じでないというのでなければ、該２つの画像ＸおよびＹは変更されずに通される。それらの画像が同じサイズでない場合、それらは、画像の一方または両方をゼロ値の画素で適切に埋込むことによって水平方向寸法および垂直方向寸法に関して同じになるように変更される（それぞれの寸法に関して、その寸法（方向）の最小の解像度で画像を埋込む）。寸法を同じにした画像は均等化基準画像Ｘ’および均等化歪画像Ｙ’と呼ぶものとする。それらはいずれもＭ×Ｎピクセルのサイズである。

以下の記載を簡潔にするために、均等化基準画像Ｘ’および均等化歪画像Ｙ’は、画素埋込が任意の所与の場合において実際には必要でない場合があるという理解の上で、引き続き、単に基準画像Ｘおよび歪画像Ｙと呼ぶことにする。画像処理の概念に詳しい人は、画素埋込が行われた場合、処理のステップ（実施形態の記載には示されない）は、レジストレーション画像Ｙ^＊から画素の行または列を削除してレジストレーション画像Ｙ^＊の寸法を（真の）基準画像Ｘの寸法に確実にマッチングさせるように行われ得ることを理解する。

Ｘ−ＳＩＦＴ抽出モジュール２４０およびＹ−ＳＩＦＴ抽出モジュール２５０において、スケール不変特徴変換（ＳＩＦＴ）の構成は、Ｄ．Ｇ．Ｌｏｗｅ「Ｄｉｓｔｉｎｃｔｉｖｅ Ｉｍａｇｅ Ｆｅａｔｕｒｅｓ ｆｒｏｍ Ｓｃａｌｅ−Ｉｎｖａｒｉａｎｔ Ｋｅｙｐｏｉｎｔｓ」、Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｊｏｕｒｎａｌ ｏｆ Ｃｏｍｐｕｔｅｒ Ｖｉｓｉｏｎ、第６０巻、第２号、第９１−１１０ページ、２００４年によって記載されるように、基準画像Ｘおよび歪画像Ｙから、それぞれキーポイントの組ＫｘおよびＫｙの形態で抽出される。

キーポイントの組ＫｘおよびＫｙはキーポイントマッチングモジュール２６０に入力される。

キーポイントマッチングモジュール２６０において、マッチングしたキーポイントサブセットＫｘ’およびＫｙ’として、所定の数のキーポイントがキーポイントの組「Ｋｘ」および「Ｋｙ」から選択され、これは、最も近いマッチング特徴ベクトルを有する。選択されたキーポイントサブセットＫｘ’およびＫｙ’は、基準画像Ｘおよび歪画像Ｙの間の幾何学的関係を最も近く近似する順アフィン変換「ＦＡＴ」のパラメータを計算するように、アフィン変換生成モジュール２７０に送られる。ＦＡＴは変換反転モジュール２８０に送られて反転され、逆アフィン変換ＩＡＴが生成される。これは、図１のアライメントパラメータの組の第１の実施形態のバージョンを表す。

アフィン変換パラメータ推定サブシステム２１０において、逆アフィン変換ＩＡＴは歪画像Ｙをレジストレーション画像Ｙ^＊に変換するために適用される。

図２ｂは、画像前処理システム１３０の第２の実施形態である、代替のアフィン変換に基づく画像アライメントシステム２００．２のブロック図である。図２ｂの代替のアフィン変換に基づく画像アライメントシステム２００．２は、変換反転モジュール２８０およびアフィン画像変換モジュール２２０が単一の、より効率的なアフィン補償画像生成モジュール２９０に結合されている点を除いて、図２ａのアフィン変換に基づく画像アライメントシステム２００．１と同じ機能を提供する。

図３は、本発明の第１の実施形態によるアフィン変換に基づく再アライメント方法３００の流れ図であり、以下のステップを含む。
３１０ 「ゼロ埋込画像」
３２０ 「重要な特徴抽出」
３３０ 「キーポイントをマッチングさせる」
３４０ 「マッチングする画素の数＜ＮＫ？」
３５０ 「アフィンパラメータを推定する」
３６０ 「変換を反転する」
３７０ 「アフィン変換を実行する」
３８０ 「アフィン補償画像を生成する」、および
３９０ 「空白の画素を埋込む」

基準画像Ｘをｘ（ｉ，ｊ）（ｉ＝１〜ｍ、ｊ＝１〜ｎ）として表し、歪画像Ｙをｙ（ｉ，ｊ）（ｉ＝１〜ｍ’、ｊ＝１〜ｎ’）として表す。基準画像Ｘは、圧縮、ノイズ汚染、ならびにあり得るアフィン変換を含む様々な歪みを受け、歪画像Ｙとなったものであり得る。２つの画像（ＸおよびＹ）が完全に異なる画像であることは考えない。

図２ａの画像埋込モジュール２３０において実行されるステップ３１０「ゼロ埋込画像」において、基準画像Ｘの寸法ｍおよびｎが、歪画像Ｙの対応する寸法ｍ’およびｎ’と比較される。

基準画像Ｘと歪画像Ｙが同一の数の画像を有するものでない場合、均等化画像ＸおよびＹの寸法をＭおよびＮに設定する。ここでＭはｍとｍ’のうち大きいほうであり、Ｎはｎとｎ’のうち大きいほうである。いずれかの画像が、オリジナル画像の周囲で寸法の小さいところにゼロ画素を追加することによって埋め込まれて、そのときの場合により、基準画像Ｘ及び歪画像Ｙを置き換える均等化基準画像Ｘ’または均等化歪画像Ｙ’を生成し得る。記載を簡潔にするために、埋込が行われたか否かにかかわらず、画像は基準画像Ｘおよび歪画像Ｙと呼ぶものとする。

図２ａのＸ−ＳＩＦＴ抽出モジュール２４０およびＹ−ＳＩＦＴ抽出モジュール２５０において実施されるステップ３２０「重要な特徴を抽出する（Extract Key Feature）」において、スケール不変特徴変換（ＳＩＦＴ）は、基準画像Ｘおよび歪画像Ｙにおいて、ローカルな（局所的な）特徴をそれぞれ検出および記述するために計算される。詳細な方法は、Ｄａｖｉｄ Ｌｏｗｅ「Ｄｉｓｔｉｎｃｔｉｖｅ Ｉｍａｇｅ Ｆｅａｔｕｒｅｓ ｆｒｏｍ Ｓｃａｌｅ−Ｉｎｖａｒｉａｎｔ Ｋｅｙｐｏｉｎｔｓ」、Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｊｏｕｒｎａｌ ｏｆ Ｃｏｍｐｕｔｅｒ Ｖｉｓｉｏｎ、第６０巻、第２号、第９１−１１０ページ、２００４年によって提案されている。

ＳＩＦＴキーポイントが基準画像Ｘおよび歪画像Ｙの両方から抽出される。Ｄａｖｉｄ Ｌｏｗｅの「ＳＩＦＴ」手法に続いて、キーポイントのリストが得られ、そこでは各キーポイントはローカルな記述子のベクトル、すなわち画像内での位置、および、色、テクスチャなどのコンテンツの視覚的特徴の記述を含む。基準画像Ｘおよび歪画像Ｙに関して、それぞれ複数のＫｘおよびＫｙキーポイント（またはその組）が、この方法によって得られる。各キーポイントは１２８の値からなる記述子ベクトルによって表され、行指標、列指標、倍率、およびキーポイントの方向という４つの値を含む位置レコードと結びつけられてもよい。

そして、キーポイント記述子は、Ｎ次元ユークリッド空間における点を表すように正規化される。ここでＮ＝１２８であり、各記述子の１２８の値に対応する。

図２ａのキーポイントマッチングモジュール２６０において実施されるステップ３３０「キーポイントをマッチングさせる」において、基準画像Ｘのキーポイントと歪画像Ｙのキーポイントの間の近接マッチング（ｃｌｏｓｅ ｍａｔｃｈ）の数ＮＫが計算される。ここでマッチングの近さは１２８次元空間の距離値によって決定される。キーポイントの特徴に関して、その距離が、２番目に近いマッチングのｄｉｓｔＲａｔｉｏ倍の距離よりも小さい場合にのみ、マッチングが認められる。パラメータｄｉｓｔＲａｔｉｏは０．２から０．６の範囲内で選択されてよい。ｄｉｓｔＲａｔｉｏとして０．４の値が以下の実験において用いられる。

計算の効率化のために、真のユークリッド距離ではなく単位ベクトル間のドット積を計算するとより効率的であり得る。単位ベクトルのドット積のａｃｏｓによって計算される角度の比は、小さい角度についてのユークリッド距離の比の良い近似であることに留意する。したがって、最も近い点のベクトル角と２番目に近い近接点のベクトル角の比がｄｉｓｔＲａｔｉｏよりも小さい場合にのみマッチングが保たれる。

ステップ３４０「マッチング画素数＜ＮＫ？」において、前のステップ３３０においてマッチングしているキーポイントの対が少なくともＮＫ個見つかったかどうかを判定する。ＮＫの値は好ましくは４に設定される。見つかったマッチング画素の対が４よりも少ない場合（ステップ３４０からＹｅｓで抜ける）、歪画像Ｙのアフィン前処理を実行することは不可能であり機能は停止しなければならない。そうでなければ（ステップ３４０からＮｏで抜ける）、実行はステップ３５０へと続く。

図２ａのアフィン変換生成モジュール２６０で実行されるステップ３５０「アフィンパラメータを推定する」において、マッチングしているキーポイントのペアの間のグローバルアフィン変換が、線形最小二乗法を用いて決定され得る。

詳細には、アフィン変換におけるパラメータは以下のように決定される。

基準画像の点［ｉ_ｋｊ_ｋ］^Ｔから歪画像の点［ｕ_ｋｖ_ｋ］^Ｔへのアフィン変換は以下のように記される。

ここでｋはマッチングしているキーポイントの指数である。未知の変数は変換パラメータａ_１、ａ_２、ａ_３、ａ_４、ｔ_１およびｔ_２である。方程式（１）は以下のように書き直される。

ここで、マッチングしているキーポイントの各組が、左端の行列の２つの行と右端のベクトルに寄与する。右端のベクトルはしたがって［ｕ_１ｖ_１ｕ_２ｖ_２．．．ｕ_ＮＫｖ_ＮＫ］^Ｔである。この方程式は線形最小二乗法で解くことができる。方程式（２）をＰＸ＝Ｑと略す。したがって、
Ｘ＝（Ｐ^ＴＰ）^−１Ｐ^ＴＱ （３）
を得る。

したがって、変換パラメータはａ_１＝Ｘ（１）、ａ_２＝Ｘ（２）、ａ_３＝Ｘ（３）、ａ_４＝Ｘ（４）、ｔ_１＝Ｘ（５）、そしてｔ_２＝Ｘ（６）である。

ステップ３５０「アフィンパラメータを推定する」の出力は、したがって、上記のようにして導かれた、アフィン変換パラメータａ１、ａ２、ａ３、ａ４、ｔ１およびｔ２の組である。それは基準画像Ｘと歪画像Ｙの間の順アフィン変換（ＦＡＴ）を示す。歪画像は、ＦＡＴの反転を歪画像に適用することにより、基準画像Ｘに似せるように、幾何学的に操作することができる。すなわち、スケーリングしたり、回転したり、シフト（平行移動）したり、剪断（ｓｈｅａｒ）したりさえもできる。それによって歪画像Ｙをレジストレーション画像Ｙ^＊、例えばアフィン補償画像ｙ^＊（ｉ^＊，ｊ^＊）に変換し、該画像は基準画像Ｘと同じ解像度を有する。

ｉ^＊∈［１，ｍ］かつｊ^＊∈［１，ｎ］として、アフィン補償画像ｙ^＊（ｉ^＊，ｊ^＊）に対して、［ｉ^＊ｊ^＊］^Ｔを２つの方法で計算することができる。

図２ａの変換反転モジュール２６０において実行されるステップ３６０「変換を反転する」において、逆アフィン変換は以下のように実施される。

その後に、図２ａのアフィン画像変換モジュール２２０において実行されるステップ３７０「アフィン変換を実行する」が行われる。それによって前処理された画像が、
ｙ^＊（ｉ^＊，ｊ^＊）＝ｙ（ｕ，ｖ） （５）
として得られる。

しかしながら、計算された値ｉ^＊およびｊ^＊は整数値ではないこともあるため、ｙ^＊（ｉ^＊，ｊ^＊）の要素、整数位置でのアフィン補償画像ｙ^＊を見つける問題に取り組む必要がある（存在するが、ステップ３７０に明示的されていない）。

図２ｂのアフィン補償画像生成モジュール２９０において実行される、好ましいステップ３８０「アフィン補償画像を生成する」において、上記のステップ３６０および３７０の機能は組み合わせられて、画素が整数指数位置に投影されて直接アフィン補償画像ｙ^＊（ｉ^＊，ｊ^＊）を生成する。画素値はバイリニア（ｂｉｌｉｎｅａｒ）補間によって以下のように計算される。

アフィン補償画像ｙ^＊（ｉ^＊，ｊ^＊）のすべての画素（ｉ^＊，ｊ^＊）に対して、以下の式が成立する。

変形された指標値（ｕ，ｖ）^Ｔが実数値であるため、バイリニア補間は値ｙ^＊（ｉ^＊，ｊ^＊）を計算するために使用可能である。

指標値ｕおよびｖが１≦ｕ≦ｍかつ１≦ｖ≦ｎの範囲にある場合、以下の方程式（７）はステップ３８０の処理の簡略表現を提供する。

ここで、

および

は、それぞれ実数値ｕの天井関数および床関数である。これに対して、ｕまたはｖがこれらの範囲外である場合、すなわちｕ＜１またはｍ＜ｕまたはｖ＜１またはｎ＜ｖである場合、対応する画素値を０に割り当てる：ｙ^＊（ｉ^＊，ｊ^＊）＝０。

ステップ３８０、すなわち方程式（６）および（７）の例示的な実施形態は、以下のＭＡＴＬＡＢコードサンプル１にさらに示される。留意点：ＭＡＴＬＡＢ（ｍａｔｒｉｘ ｌａｂｏｒａｔｏｒｙ）はＭａｔｈ−Ｗｏｒｋｓ（ｈｔｔｐ：／／ｗｗｗ．ｍａｔｈｗｏｒｋｓ．ｃｏｍ）によって開発された数値計算環境および第４世代プログラミング言語である。ＭＡＴＬＡＢコードサンプル１。

ステップ３９０「任意的に空白の画素を埋める」において、アフィン補償画像ｙ^＊の画像境界近くの任意の空白の画素は空白のままにされてもよく、任意的に、画像境界の非空白画像画素のミラーリング値によって埋込まれてもよい。

例えば、以下の画像がある場合、

この場合、欠けた画素（ゼロ値）のミラーリング値による埋込みは、空白値と非空白値の間の境界における非空白画素値をミラーリングすることによって、実行される。これは以下を与える。

本発明の第１の実施形態の評価の目的のため、これらの画素は空白のままとする。

要するに、アフィン変換に基づく再アライメント方法３００は、基準画像Ｘと歪画像Ｙを入力として受信し、アフィン変換技術と結合したＳＩＦＴ技術を用いて、アライメントされたレジストレーション画像Ｙ^＊を出力する。

図４ａ、図４ｂおよび図４ｃは図３のアフィン変換に基づく再アライメント方法３００の画像的な結果を示す。基準画像Ｘの代わりに用いられるオリジナルの画像が図４ａに示される。アフィン変換に基づく再アライメント方法３００を説明するために、図４ａのオリジナル画像のアフィン変換によって、歪画像Ｙとして使用するための１つの歪画像が作成され、スケーリング、回転、平行移動、および剪断の歪みの組み合わせを呈する。歪画像Ｙは図４ｂに示される。図４ｂの歪画像Ｙは、図３のアフィン変換に基づく再アライメント方法３００により前処理され、前処理画像、すなわち図４ｃに示すレジストレーション画像Ｙ^＊を生成した。視覚的な検査によって、前処理方法は歪画像の基準画像へのアライメントにおいて優れた仕事を行うことが可能であることがわかる。

アフィン変換に基づく再アライメント方法３００の性能をテストするために、Ｈ．Ｒ．Ｓｈｅｉｋｈ、Ｚ．Ｗａｎｇ、Ｌ．ＣｏｒｍａｃｋおよびＡ．Ｃ．Ｂｏｖｉｋにより刊行され、ｈｔｔｐ：／／ｌｉｖｅ．ｅｃｅ．ｕｔｅｘａｓ．ｅｄｕ／ｒｅｓｅａｒｃｈ／ｑｕａｌｉｔｙからインターネット上で利用可能な、ＬＩＶＥ画質評価データベースリリース２からの画像を用いて多くの実験が実施された。このデータベースは、５つのタイプの歪みを用いた２９のオリジナルの画像に由来する７７９の歪画像からなる。歪みは、ＪＰＥＧ圧縮、ＪＰＥＧ２０００圧縮、ガウスホワイトノイズ（ＧＷＮ）、ガウスぼかし（ＧＢｌｕｒ）、およびレイリー高速フェーディング（ＦＦ）チャネルモデルを含む。加えて、歪画像は、スケーリング、回転、および水平垂直両方の方向の空間的シフトによってオリジナルの画像から生成された。

テストのために、オリジナル（基準）の画像にアライメントされた、前処理されたレジストレーション画像Ｙ^＊が、図１（画像評価プロセス１６０参照）に示されるように、視覚品質評価をされた。基準画像Ｘと前処理されたバージョンの歪画像Ｙ、すなわちレジストレーション画像Ｙ^＊の視覚的質の測定には、標準的な測定基準が用いられた。３つの完全参照測定基準（ｆｕｌｌ−ｒｅｆｅｒｅｎｃｅ ｍｅｔｒｉｃｓ）であるピーク信号ノイズ比（ＰＳＮＲ）、平均構造的類似性指標測定（ＭＳＳＩＭ）、および視覚情報忠実度（ＶＩＦ）測定もまた、画像の視覚的質の測定のために考慮された。

３つの性能測定法が考慮された。すなわち、差分平均オピニオン評点（ＤＭＯＳ）と非線形回帰後の客観的モデル出力の間の相関係数（ＣＣ）、二乗平均平方根誤差（ＲＭＳＥ）、スピアマン順位相関係数（ＳＲＯＣＣ）である。ＳＲＯＣＣは以下のように定義される。

ここでＮは評価が行われる画像内の画素の数であり、ｄ_ｉは主観的評価および客観的評価におけるｉ番目の画像の順位の差である。Ｚ．ＷａｎｇおよびＱ．Ｌｉ「Ｉｎｆｏｒｍａｔｉｏｎ ｃｏｎｔｅｎｔ ｗｅｉｇｈｔｉｎｇ ｆｏｒ ｐｅｒｃｅｐｔｕａｌ ｉｍａｇｅ ｑｕａｌｉｔｙ ａｓｓｅｓｓｍｅｎｔ」、ＩＥＥＥ Ｔｒａｎｓａｃｔｉｏｎｓ ｏｎ Ｉｍａｇｅ Ｐｒｏｃｅｓｓｉｎｇ、第２０巻、第５号、第１１８５−１１９８ページ、２０１１年に記載のように、ＳＲＯＣＣはノンパラメトリックの順位に基づく相関測定基準であり、主観的評点と客観的評点の間のいかなる単調な非線形マッピングとも独立している。

以下の表１−表５は、前処理なし（「オリジナル」方法とする）の歪画像の視覚的質測定と、本発明の第１の実施形態によるアフィン前処理（「提案」方法とする）による歪画像の視覚的質測定との比較から、実験結果をまとめたものである。

表１は３つの完全参照測定基準、すなわちＰＳＮＲ、ＭＳＳＩＭおよびＶＩＦについていかなる歪みも追加されない、このデータベースの実験結果を示す。最良の実施方法は太いフォントで強調する。表１に示された結果は、歪画像の前処理が、オリジナルのデータベースに関するオリジナルの方法よりも画質評価の特性に悪影響を与えないことを示す。すなわち、実際により良くなっている。

表２は、ＬＩＶＥデータベースのすべての画像が、倍率０．９のスケーリング、回転角０．１×１８０／π度の回転、水平方向に４画素、および垂直方向に４画素の空間シフトにより歪んだ場合の結果を表に示す。この表において、標準の完全参照測定基準、すなわちＰＳＮＲ、ＭＳＳＩＭおよびＶＩＦが非常に低い質評点を得る。それに対して、前処理方法は非常に良好な視覚的質評点を与える、すなわち、前処理された画像について行われた客観的質評価は、前処理されない画像よりもはるかに主観的評価との相関が高い。最もうまくいく方法は太いフォントで強調する。実際に、前処理なしでは、既存の質測定基準は非常にパフォーマンスが低いことがわかる。

表２は、歪みがスケーリング、回転およびシフトを含む場合、オリジナルの方法に対してパフォーマンスが向上することを示す。そのようなパフォーマンスの向上は、図８、図９および図１０に示された、提案の測定基準と人の評点（ｈｕｍａｎ ｓｃｏｒｅ）の間の相関がより良いことから見てとれる。図１１、図１２および図１３に示された結果は、本発明の提案の方法はオリジナルの方法よりも、様々な歪みに対してセンシティブではないことを示す。

歪画像の提案の前処理に関する目的および根拠が以下に示される。視覚的質にはまったく影響を与えるべきでない、例えば回転、スケーリングなどのアフィン変換を補償することが求められている。事実、図１１および図１２を見ると、提案の方法は小さな回転およびシフトにはかなりフラットである（不変である）ことがわかる。しかし、画像がかなりスケーリングされて、それによって画像の細部が失われる場合、質を同じに保つことが期待できないため、スケーリングはより扱いにくくなる。しかし、ＶＩＦに関して図１３ｃに示すように、この方法は、倍率の変化で緩やかに減少する合理的な評点を与えることが望まれる。この場合には、提案の方法の評点は緩やかに減少して合理的な状態にとどまる（例えばスケーリングが０．５の場合に評点０．７）。これに対してオリジナルの方法（すなわち前処理なし）は、スケーリングが０．５のとき評点が０．１未満となり、オリジナルの解像度の半分に画像をスケーリングした場合に質がそこまで低くはならないであろうため、直感的にあまり意味がない。

表１．ＬＩＶＥ画像データベース内のすべての歪テスト画像に関する画質評価の全体的なパフォーマンス比較。

表２．ＬＩＶＥ画像データベース内の、スケーリング、回転、および水平方向及び垂直方向両方の空間的シフトによって歪みを受けた、各画像に関する画質評価の全体的なパフォーマンス比較。

表３−表５はＰＳＮＲ、ＭＳＳＩＭおよびＶＩＦとＬＩＶＥ画像データベースの５種類の歪みとの対応に関してＣＣ、ＳＲＯＣＣおよびＲＭＳＥの値をそれぞれ示す。アフィン前変換があるときの評点はそれがないときの評点よりも良い。ＶＩＦに関しては、アフィン前処理があるときの評点はそれがないときの評点と同程度である。
表３．ＬＩＶＥ画像データベースに関する非線形回帰後のＣＣ値。

表４。ＬＩＶＥ画像データベースに関する非線形回帰後のＳＲＯＣＣ値。

表５。ＬＩＶＥ画像データベースに関する非線形回帰後のＲＭＳＥ値。

図５から図１３はさらに詳細に実験結果を描写する結果の画像を示す。

図５ａ、図６ａおよび図７ａは、アフィン歪みなしの「オリジナル」の測定（すなわち、アフィン変換パラメータ推定およびアフィン補償がない）に関する、ＬＩＶＥ画像データベースのすべての歪画像に対する、ＤＭＯＳ（縦軸）と、それぞれＰＳＮＲ、ＭＳＳＩＭおよびＶＩＦの評点（横軸）の散布図である。これに対して図５ｂ、図６ｂおよび図７ｂは、画像のアフィン前処理後の測定に関する、対応する散布図である。

図８ａ、図９ａおよび図１０ａは、追加のアフィン歪み（倍率０．９のスケーリング、回転角０．１×１８０／π度の回転、水平方向に４画素、および垂直方向に４画素の空間シフト）が加えられた、「オリジナル」の測定に関する、ＬＩＶＥ画像データベースのすべての歪画像に対する、ＤＭＯＳ（縦軸）と、それぞれＰＳＮＲ、ＭＳＳＩＭおよびＶＩＦの評点（横軸）の散布図である。これに対して図８ｂ、図９ｂおよび図１０ｂは、画像のアフィン前処理後の測定に関する、対応する散布図である。

図８から図１０の散布図は、追加の歪みが加えられた画像の相関を示し、オリジナルの測定と比較しての前処理による顕著な改善を示す。これは水平軸のスケーリングを異なるものにすることにより一層明らかになる。前処理なしの画像の場合、ＶＩＦ評点は０．００４から０．０２のような低さである。しかしながら、アフィン前処理後は、ＶＩＦ評点は０と１の間の範囲にある。

図１１ａ、図１１ｂ、および図１１ｃは、それぞれ測定基準ＰＳＮＲ、ＭＳＳＩＭおよびＶＩＦに関する、ＬＩＶＥデータベースの選択された画像を平行移動によって歪ませた実験の結果を示すグラフである。平行移動は１から１５画素の範囲にある。

同様に、図１２ａ、図１２ｂ、および図１２ｃは、それぞれ測定基準ＰＳＮＲ、ＭＳＳＩＭおよびＶＩＦに関する、ＬＩＶＥデータベースの選択された画像を１度から１５度までの回転によって歪ませた実験の結果を示すグラフである。図１３ａ、図１３ｂ、および図１３ｃは、それぞれ測定基準ＰＳＮＲ、ＭＳＳＩＭおよびＶＩＦに関する、ＬＩＶＥデータベースの選択された画像を０．５から１．０までの倍率によって歪ませた実験の結果を示すグラフである。

図１１ａ、図１１ｂ、図１２ａ、および図１２ｂにおける水平軸のスケールは、比較をより明確に示すために（０，０）を含まない。実際、（０，０）においては歪みはなく、したがってオリジナルの画像およびアフィン前処理された画像の両方が同じ測定基準の評点となるであろう。

平行移動歪みに関するグラフに見られるように、歪画像をアフィン前処理する提案の方法は、測定基準ＰＳＮＲ、ＭＳＳＩＭ、およびＶＩＦに関してほぼ一定の評点を出す。小さな平行移動は視覚的質に影響を与えるべきでないため、ほぼ一定の評点は、良好な前処理が実行される場合に期待される。回転歪みに関しては、提案の方法はＭＳＳＩＭおよびＶＩＦに関してほぼ一定の評点を出す。ＰＳＮＲ値は空間シフトが増加するにつれて減少するが、提案の方法は前処理なしの標準的なＰＳＮＲよりはなお良好である。また、小さな回転は視覚的質に影響を与えるべきでないため、ほぼ一定の評点は、良好な前処理が実行される場合に期待される。測定基準の評点が平行移動および回転においてほとんど一定のままであるという事実は、この測定基準を適用する前の提案の前処理の有効性を示す。スケーリングの歪みに関して、提案の方法はＰＳＮＲ、ＭＳＳＩＭ、およびＶＩＦに関する標準の測定基準よりも顕著に高い評点を出す。倍率が減少するにつれて質の低下が予想されるが、これは、画像の品質は、小さくなるにつれて、低下する（一層ぼやける）ためである。しかしながら、オリジナルの方法において示されるように激しく減少することは直感的にありえない。提案の方法において示されるように、質の緩やかな低下がより適切である。

これは図７ａ、図７ｂおよび図７ｃに示されている。そこでは測定基準による値は倍率が減少するにつれて減少する。しかし変化は本発明の前処理方法を用いた場合は緩やかであり、それは、画像のスケーリングに関連した質の緩やかな減少を表す。出願人が発見したことは、画像がアフィン歪みを受ける場合にはＭＳＳＩＭおよびＶＩＦはＰＳＮＲよりも良好な測定基準であるということである。しかしながら、本発明の前処理ステップで生成された３つの測定基準すべてが、前処理なしのこれらの標準的な測定基準よりも高い評点を出す。提案の前処理方法は異なる種類の歪みに対して非常に堅牢であることが容易にわかる。

Ｍ．Ｐ，Ｓａｍｐａｔ、Ｚ．Ｗａｎｇ、Ｓ．Ｇｕｐｔａ、Ａ．Ｃ．ＢｏｖｉｋおよびＭ．Ｋ．Ｍａｒｋｅｙらの「Ｃｏｍｐｌｅｘ ｗａｖｅｌｅｔ ｓｔｒｕｃｔｕｒｅ ｓｉｍｉｌａｒｉｔｙ：Ａ ｎｅｗ ｉｍａｇｅ ｓｉｍｉｌａｒｉｔｙ ｉｎｄｅｘ」、ＩＥＥＥ Ｔｒａｎｓａｃｔｉｏｎｓ ｏｎ Ｉｍａｇｅ Ｐｒｏｃｅｓｓｉｎｇ、第１８巻、第１１号、第２３８５−２４０１ページ、２００９年においてＳａｍｐａｔ等は、小さい平行移動および小さい回転にのみ不変な、複素ウェーブレット領域における新しい画像類似度指標（ＣＷ−ＳＳＩＭ）を提示した。これに対し、本特許出願において提案の前処理方法は、任意量の平行移動、回転、および倍率に対してより堅牢である。このことは、本発明の提案の方法が、異なる種類の歪みのもとでの画像の視覚的な質を測定するのに適していることを示す。

画像前処理システム１３０の第１および第２の実施形態が、ＳＩＦＴキーポイントマッチングおよびアフィン変換を用いた「アフィン変換に基づく」実施形態であるのに対し、ラドン変換に基づく追加の実施形態が以下に記載される。

アフィン変換に基づく方法とラドン変換に基づく方法のいずれも、歪画像を基準画像にアライメントする処理を行うという同じ目的を果たすが、能力にいくつかの相違がある。例えば、アフィン変換に基づく方法はスキューおよびフリップ（反転）を含んだ多くの幾何学的歪みを補正することができるが、ラドン変換に基づく方法はスキューおよびフリップの歪みを補正することができない。その一方で、非常にノイズが多い画像が存在する場合のラドン変換に基づく方法の堅牢性は、アフィン変換に基づく方法の性能を上回る。

画像前処理システム１３０のラドン変換に基づく実施形態は、歪画像Ｙを基準画像Ｘによって「レジストレーションする」ための倍率および回転角を決定する方法に基づく。これらの追加実施形態は効率的な実施につながり、ノイズの存在時に優れて堅牢であるが、上記第１および第２の実施形態のほうが優れているある種の歪み（例えば剪断歪み。図４ｂ参照）に関して制限される。これらの実施形態はさらに、同じ倍率が水平方向と垂直方向の両方にあてはまるという前提に基づいている。しかしながら、提案の方法はノイズが多い環境では非常に良好にはたらく。

図１４は、以下の機能を実行する複数のブロックを含むラドン変換に基づくアライメントシステムの構成ブロックの組のダイアグラム１４００を示す。

入力画像、すなわち基準画像Ｘと歪画像Ｙの水平方向寸法を均等化する、寸法均等化機能１４１０。

画像をその重心にセンタリングする、センタリング機能１４２０。

歪画像Ｙのサイズに対する基準画像Ｘのサイズの比として倍率「ａ」を決定する、倍率推定機能１４３０。

基準画像Ｘと歪画像Ｙの方向の間に存在する回転角θ_０を決定する、角度推定機能１４４０。

歪画像Ｙまたはそれを倍率「ａ」によって処理したバージョンをスケーリングして、そのサイズを基準画像Ｘのサイズにほぼ等しくする、スケーリング機能１４５０。

基準画像Ｘまたはそれを倍率「ａ」の逆数によって処理したバージョンをプレスケーリングして、そのサイズを歪画像Ｙまたはそれを処理したバージョンのサイズにほぼ等しくする、前スケーリング機能１４６０。

歪画像Ｙまたはそれの処理したバージョンを回転角θ_０だけ回転する、回転機能１４７０。

歪画像Ｙの処理したバージョンの平行移動を行い、基準画像Ｘまたはそれの処理したバージョンに垂直方向および水平方向でアライメントする、平行移動機能１４８０。

以下に記載される本発明の実施形態は、歪画像Ｙと基準画像Ｘを処理することによってレジストレーション画像Ｙ^＊を生成するように、様々な組み合わせおよび順番で構成要素の組１４００を使用する。

図１５は図１の画像前処理システム１３０の第３実施形態である、第１のラドン変換に基づく画像アライメントシステム１５００のブロック図である。それは以下のモジュールを含む。

倍率「ａ」を決定する、倍率推定モジュール１５０２。

歪画像Ｙからマスキングされた画像Ｙ_０を抽出する、画像マスキングモジュール１５０４。

マスキングされた画像Ｙ_０からセンタリングパラメータ「ｃ_ｒ」および「ｃ_ｃ」を決定する、重心計算モジュール１５０６。

センタリングパラメータ「ｃ_ｒ」および「ｃ_ｃ」を用いて歪画像Ｙをセンタリングして、センタリングされた画像Ｙ_１を生成する、画像センタリングモジュール１５０８。

倍率「ａ」を用いてセンタリングされた画像Ｙ_１をサイズ変更して、サイズ変更された歪画像Ｙ_２とする、画像サイズ変更モジュール１５１０。

サイズ変更された歪画像Ｙ_２は基準画像Ｘと同じサイズであり、センタリングされているが、さらに回転する必要がある場合がある。

第１のラドン変換に基づく画像アライメントシステム１５００はさらに以下を含む。

サイズ変更された歪画像Ｙ_２と基準画像Ｘのラドン変換Ｒ_１およびＲ_２をそれぞれ生成するラドン変換モジュール１５１２、およびラドン変換Ｒ２とＲ１から回転角「θ_０」を推定するラドン変換相関モジュール１５１６を含む、回転角決定ユニット１５１１。

サイズ変更された歪画像Ｙ_２が回転角「θ_０」だけ回転され、補正された画像Ｙ_３を形成するが、当該画像は基準画像Ｘからさらに横方向にオフセットしていてもよい、画像回転モジュール１５１８。

補正された画像Ｙ_３と基準画像Ｘの間にオフセットベクトル「ＴＶ」が決定される、平行移動推定モジュール１５２０。

オフセットベクトル「ＴＶ」が補正された画像Ｙ_３に適用されてレジストレーション画像Ｙ^＊を生成する、画像平行移動モジュール１５２２。

第３の実施形態の方法に関して、歪画像Ｙ内の、実際の視覚的な内容（コンテンツ）の周囲に黒い背景があることが前提となっている。回転角が０度、９０度、１８０度および２７０度から顕著に離れるにつれ、多くの画素がこの方法においては失われることに留意することも重要である。それらの角度の値から離れるほど、失われる画素の数が多くなる。したがって、アライメント方法の後に画像評価方法が続く場合、これらの角度の値から離れるほど、質の評価はより不正確になる。

図１６は、第１のラドン変換に基づく画像アライメントシステム１５００のモジュールにおいて例示されることができる、第１のラドン変換に基づく方法１６００を示し、それは以下のステップを含む。
３１０ 「ゼロ埋込画像」
１６１０ 「倍率ａを推定する」
１６２０ 「画像Ｙの重心を計算する」
１６３０ 「画像ＹをセンタリングしてＹ_１を生成する」
１６４０ 「画像Ｙ_１をサイズ変更してＹ_３を生成する」
１６５０ 「ラドン変換Ｒ１およびＲ２を実行する」
１６６０ 「変換を相関してθ_０を生成する」
１６７０ 「画像Ｙ_１を回転してＹ_３を生成する」および
１６８０ 「補正された画像Ｙ_３を平行移動する」

前記のステップ３１０「ゼロ埋込画像」で基準画像Ｘと歪画像Ｙのサイズを均等化した後、以下の一連のステップが行われ、そこでは、歪画像Ｙがスケーリングされ、回転され、平行移動され、レジストレーション画像Ｙ^＊となる。

倍率推定モジュール１５０２において実行される、ステップ１６１０「倍率ａを推定する」において、フーリエ変換を用いて行うことができる各画像の輝度（ｌｕｍｉｎａｎｃｅ）の平均化と、その結果を割ることによって、倍率「ａ」が決定される。これは、歪画像の生成において、平均の輝度を変更するいかなる操作（例えばコントラストまたはブライトネスの操作など）も行われないことを前提とする。

歪画像Ｙを、基準画像Ｘを回転してスケーリングしたものとする。
それぞれサイズが（Ｍ×Ｎ）である、画像ＸおよびＹのフーリエ変換Ｆ１およびＦ２は、以下の方程式（９）及び（１０）によって与えられる。

極座標系において、これらの大きさは方程式（１１）において与えられるように関連する。Ｂ．Ｓ．ＲｅｄｄｙおよびＢ．Ｎ．Ｃｈａｔｔｅｒｊｉ「Ａｎ ＦＦＴ−Ｂａｓｅｄ ｔｅｃｈｎｉｑｕｅ ｆｏｒ ｔｒａｎｓｌａｔｉｏｎ，ｒｏｔａｔｉｏｎ ａｎｄ ｓｃａｌｅ−ｉｎｖａｒｉａｎｔ ｉｍａｇｅ ｒｅｇｉｓｔｒａｔｉｏｎ」、ＩＥＥＥ Ｔｒａｎｓａｃｔｉｏｎｓ ｏｎ Ｉｍａｇｅ Ｐｒｏｃｅｓｓｉｎｇ、第５巻、第８号、第１２６６−１２７１ページ、１９９６年参照。詳細は、

である。
ここでａは倍率（水平方向および垂直方向に同じ倍率を適用する）であり、θ_０は２つの画像の間の回転角であり、

θ＝ｔａｎ^−１（ｎ／ｍ） （１３）
である。
したがって、倍率は下記のように得られる。

倍率ａは画像ＸおよびＹの２Ｄ ＦＦＴを計算することなく、より簡潔に得ることもできる。なぜなら

および

だからである。

ここで、歪画像Ｙがスケーリング同様に平行移動および回転されてもよいと考える。シフト操作が歪画像Ｙの画素の合計数に対して非常にわずかな数の画素しか取り除かない限り、倍率ａを計算するために使用された手法はなおも有効である。

画像マスキングモジュール１５０４において実行されるステップ１６２０「マスクを作成する」において、マスク画像Ｙ_０が以下の規則によって作成される。Ｙ（ｍ，ｎ）＞τの場合にはＹ_０（ｍ，ｎ）＝１、そうでなければＹ_０（ｍ，ｎ）＝０。本発明のこの実施形態において、発明者はτ＝４０を選んだ。

ステップ１６３０「画像をセンタリングする」において、画像Ｙはセンタリングされ、センタリングされた画像Ｙ_１となる。マスク画像Ｙ_０の重心（ｃ_ｒ、ｃ_ｃ）は以下のように計算される。

ここでパラメータｃ_ｒおよびｃ_ｃは、「重心」とされる画素の、それぞれ行および列の座標を表す。この演算は重心計算モジュール１５０６で行われる。

Ｙ（ｃ_ｒ，ｃ_ｃ）に対応する位置は画像Ｙの中央に移動され、平行移動は画像センタリングモジュール１５０８において行われる。センタリングされた画像Ｙ_１は方程式１９によって計算される。

ここでｍ∈［１，Ｍ］かつｎ∈［１，Ｎ］。（ｍ−ｃ_ｒ＋Ｍ／２，ｎ−ｃ_ｃ＋Ｎ／２）が指標範囲［１，Ｍ］×［１，Ｎ］の外にある場合、Ｙ_１（ｍ，ｎ）は０に設定される。

画像サイズ変更モジュール１５１０において実行される、ステップ１６４０「画像をスケーリングする」において、倍率「ａ」は、センタリングされた画像Ｙ_１を正規化された画像Ｙ_２にサイズ変更して、基準画像Ｘと同じスケールとなるように使用される。これは回転角を推定する準備において行われる。正規化された画像Ｙ_２は、Ｙ_２（ｍ，ｎ）＝Ｙ_１（ｍ／ａ，ｎ／ａ）と定義される。

正規化された画像Ｙ_２の指標ｍおよびｎは整数値でなければならないが、サイズ変更の計算の後は、それらは一般に浮動小数点値であり、切り捨てるか四捨五入して整数値にする必要がある。サイズ変更された画像内の画素値に関して、最も近い画素値を用いてもよく、双一次補間法（ｂｉｌｉｎｅａｒ ｉｎｔｅｒｐｏｌａｔｉｏｎ）を用いて近接の画素の数から画素値を決定してもよい。

基準画像Ｘと正規化された画像Ｙ２の間の回転角θ_０は、以下の２つのステップで対応するラドン変換から得られる。

ステップ１６５０「ラドン変換を実行する」は、ラドン変換モジュール１５１２において実行され、基準画像Ｘの対応するラドン変換および正規化された画像Ｙ_２の対応するラドン変換を計算する。

基準画像Ｘのラドン変換Ｒ_１は以下のように定義される。

そして、正規化された画像Ｙ_２のラドン変換Ｒ_２は以下のように定義される。

ここでδ（）はＤｉｒａｃのデルタ関数である。

Ａ．Ａｖｅｒｂｕｃｈ、Ｒ．Ｒ．Ｃｏｉｆｍａｎ、Ｄ．Ｌ．Ｄｏｎｏｈｏ、Ｍ．Ｉｓｒａｅｌｉ、およびＹ．Ｓｈｋｏｌｎｉｓｋｙ「Ａ ｆｒａｍｅｗｏｒｋ ｆｏｒ ｄｉｓｃｒｅｔｅ ｉｎｔｅｇｒａｌ ｔｒａｎｓｆｏｒｍａｔｉｏｎｓ Ｉ − ｔｈｅ ｐｓｅｕｄｏ−ｐｏｌａｒ Ｆｏｕｒｉｅｒ ｔｒａｎｓｆｏｒｍ」、ＳＩＡＭ Ｊｏｕｒｎａｌ ｏｎ Ｓｃｉｅｎｔｉｆｉｃ Ｃｏｍｐｕｔｉｎｇ、３０（２）、第７６４−７８４ページ、２００８において指摘されているように、現代の応用においては離散２Ｄラドン変換を有することが重要であり、離散２Ｄラドン変換は最近２０年間で多くの著者の注目の対象となっている。最近まで、ラドン変換は２Ｄ離散画像のコヒーレントな離散鮮明度を欠いていた、それは代数的に正確で、可逆で、すばやく計算可能である。それゆえ、Ａｖｅｒｂｕｃｈ ｅｔ ａｌ．は離散２Ｄ画像のための２Ｄ離散ラドン変換の観念を定義した。それは絶対値が１未満の勾配の線に沿った合算に基づいている。格子点でない位置における値はゼロ埋込された格子点上の三角補間を用いて定義される。離散ステップがゼロに近づくにつれて、それが連続ラドン変換に収束することから、彼らは、これらの定義が連続体の忠実な記述を提供することを証明した。離散ラドン変換に関するより詳細な記載はＡｖｅｒｂｕｃｈ ｅｔ ａｌ．にある。Ａｖｅｒｂｕｃｈ ｅｔ ａｌ．に従い、出願人はそれぞれＲ_１（ｒ，θ）およびＲ_２（ｒ，θ）と表される、ＸおよびＹ２の離散ラドン変換を行うことができる。それらは共に寸法がＫ×Ｌである。

２つの画像（基準画像Ｘおよび正規化された歪画像Ｙ_２）が同一であるか、それらの間の回転およびありうるノイズを除いて少なくとも非常に似ているという前提に基づけば、それらのラドン変換Ｒ_１（ｒ，θ）およびＲ_２（ｒ，θ）は、実質的に同じ、または、θの方向に沿っての円周方向の移動を除いて正確に同じ、すなわちＲ_２（ｒ，θ）＝Ｒ_１（ｒ，θ＋θ_０）、であるということにつながる。Ｒ_１およびＲ_２のサイズは、半径「ｒ」と角ｔｈｅｔａに関して所望の精度に基づいて選択されるパラメータＫおよびＬによって特徴づけられる。ラドン変換に基づくシステムの実施形態において、パラメータＫは用いられるソフトウェア（上記のＭＡＴＬＡＢ）の「ラドン」関数からデフォルトとして選択される。それはおよそＫ＝ｓｑｒｔ（Ｍ＊Ｍ＋Ｎ＊Ｎ）＋５、すなわち、画素数において画像の対角線の長さをわずかに上回る画像サイズに関連する。角度の精度に関して、Ｌの範囲はＬ＝０から１７９の範囲で選択され、これは約１度の精度を提供する。

ラドン変換関数は以下のようにＭＡＴＬＡＢの文書に記載される。

「Ｒ＝ｒａｄｏｎ（Ｉ、ｔｈｅｔａ）は、角度ｔｈｅｔａに対して画像強度（ｉｎｔｅｎｓｉｔｙ ｉｍａｇｅ）Ｉのラドン変換Ｒを返す。ラドン変換は特定の角度に向けた放射線に沿った画像強度の射影である。ｔｈｅｔａがスカラーである場合、Ｒはｔｈｅｔａ度に対するラドン変換を含む列ベクトルである。ｔｈｅｔａがベクトルである場合、Ｒは、各列がｔｈｅｔａ内の角度の１つに対するラドン変換である行列である。」

理論によれば、ある画像が別の画像を回転したものである場合、該２つの画像は、角度次元に沿った円周方向シフトを除いて、同じラドン変換を有することとなる。今回の場合は、基準画像Ｘと正規化された歪画像Ｙ_２はより多くの相違点を有する場合があり、変換の同一性は完全には成立しない。この理由のため、いくつかの値ｒに関して角度を計算し、角度の中央値をとることが、角度の推定をより堅牢なものとする。

ラドン変換相関化モジュール１５１６において実行される、ステップ１６６０「変換を相関させる」において、２つのラドン変換が相関方法によって体系的に比較され、歪画像Ｙの基準画像Ｘに対する最も確からしい傾き角度を推定する。

θ_０を計算するために円周方向相互相関を使用し得るが、これは時間を消費する。Ｒ_１およびＲ_２のすべての行に関して、円周方向シフトＲ_２を行い、その後に該２つの行の相互相関を計算する必要があるため、Ｏ（ＫＬ^２）の計算複雑性がある。これはＯ（Ｌ^２）の計算複雑性を有する。全部でＫ個の行があるため、全部の計算複雑性はＯ（ＫＬ^２）である。

より効率的な代替例として、１次元（１Ｄ）ＦＦＴを代わりに用いることもできる。それはＯ（ＫＬｌｏｇＬ）の計算複雑性を有する。

以下の段落はＦＦＴによる相互相関を計算する一般的な方法を提供する。そこには第３の実施形態の好ましい方法がある。

出願人は相互相関およびその高速実施の定義をここで簡潔に述べる。２つの離散した実数値関数ｆ［ｎ］およびｇ［ｎ］の間の相互相関は、Ｄ．Ｎｉｋｏｌｉｃ、Ｒ．Ｃ．Ｍｕｒｅｓａｎ、Ｗ．ＦｅｎｇおよびＷ．Ｓｉｎｇｅｒ「Ｓｃａｌｅｄ ｃｏｒｒｅｌａｔｉｏｎ ａｎａｌｙｓｉｓ：ａ ｂｅｔｔｅｒ ｗａｙ ｔｏ ｃｏｍｐｕｔｅ ａ ｃｒｏｓｓ−ｃｏｒｒｅｌｏｇｒａｍ」、Ｅｕｒｏｐｅａｎ Ｊｏｕｒｎａｌ ｏｆ Ｎｅｕｒｏｓｃｉｅｎｃｅ、第１−２１ページ、２０１２年によれば、

である。
ｆ［ｎ］とｇ［ｎ］の畳み込みはＶ．Ｓｏｂｏｌｅｖ「Ｃｏｎｖｏｌｕｔｉｏｎ ｏｆ ｆｕｎｃｔｉｏｎｓ」、Ｍｉｃｈｉｅｌ Ｈａｚｅｗｉｎｋｅｌ、Ｅｎｃｙｃｌｏｐｅｄｉａ ｏｆ Ｍａｔｈｅｍａｔｉｃｓ、Ｓｐｒｉｎｇｅｒ、ＩＳＢＮ９７８−１−５５６０８−０１０−４において以下のように与えられている。

関数ｆ［ｎ］およびｇ［ｎ］の相互相関はｈ［ｎ］＝ｆ［−ｎ］とｇ［ｎ］の畳み込みと同じである。

畳み込み理論に似て、相互相関は以下を満たす。

ここでＦＦＴは高速フーリエ変換、ｃｏｎｊは複素共役、そしてドットは成分ごとの乗法を意味する。逆ＦＦＴをとることによって、ｆ＊ｇの高速実施が可能となる。

図１７は、ステップ１６６０「変換を相関させる」の展開を示し、それは以下のステップを含む。
１７１０「Ｒ１およびＲ２の１Ｄ順ＦＦＴを計算してＲ３を生成する」
１７２０「Ｒ３／｜Ｒ３｜に１Ｄ逆ＦＦＴを適用してＲ４を生成する」
１７３０「行に関して最も大きい相互相関係数を決定する」
１７４０「配列指標最大相互相関係数の中央値からθ_０を推定する」

Ｒ_１およびＲ_２のｍ番目の行をそれぞれｒ_１（ｎ）およびｒ_２（ｎ）とし、ここでｎは、Ｒ１およびＲ２の円周方向シフトである角度の数、１からＬである。

ステップ１７１０「Ｒ１およびＲ２の１Ｄ順ＦＦＴを計算してＲ３を生成する」において、１Ｄ順ＦＦＴはｒ_１（ｎ）およびｒ_２（ｎ）のそれぞれに関して計算され、第１の中間結果Ｒ３＝ｒ_３（ｕ）を生成する。

ここで「ｃｏｎｊ」は複素共役を表し、・は成分ごとの積を表す。

ステップ１７２０「Ｒ３／｜Ｒ３｜に１Ｄ逆ＦＦＴを適用してＲ４を生成する」において、１Ｄ逆ＦＦＴはｒ_３（ｕ）／｜ｒ_３（ｕ）｜に適用される。すなわち、第２の中間結果Ｒ４＝ｒ４（ｎ）が計算される。

ｒ_２（ｎ）＝ｒ_１（ｎ＋θ_０）であるから、

ｒ_４（ｎ）＝δ（ｎ＋θ_０） （２４）
したがって、ｒ_４（ｎ）は−θ_０でピーク値１．０をとり、その他のすべての値で０をとる。Ｒ４＝ｒ４（ｎ）の値のそれぞれがＲ１とＲ２の間の、２つのラドン変換の間の与えられたオフセット角に関する相互相関係数のベクトル（すなわち配列）を表す。

ステップ１７３０「行に関して最も大きい相互相関係数を決定する」において、Ｒ４＝ｒ_４（ｎ）のすべての行から最大値の位置指標θ^＊（ｍ）を見つけ出す。すなわち、

である。

ステップ１７４０「最大相互相関係数の配列指標の中央値からθ_０を推定する」において、これらの最大値を示すそれぞれの配列指標の中央値をとることにより、２つの画像（基準画像Ｘとセンタリングされた画像Ｙ_１）の間の回転角の推定をθ_０＝ｍｅｄｉａｎ（θ^＊（ｍ））として得る。

以下のＭＡＴＬＡＢコードサンプル２が回転角の推定を実施する。

ここで図１６の記載に戻る。
画像回転モジュール１５１８において実施されるステップ１６７０「画像Ｙ_２を回転する」において、サイズ変更された画像Ｙ_２は方向差を補償するように−θ_０度だけ回転され、そして補正された画像Ｙ_３を得る。

画像平行移動モジュール１５２２において実施されるステップ１６８０「補正された画像Ｙ_３を平行移動する」において、回転された歪画像Ｙ_３は最終位置に平行移動され、それによりレジストレーション画像Ｙ^＊を生成する。

平行移動のために、Ｇ．ＶａｒｇｈｅｓｅおよびＺ．Ｗａｎｇ「Ｖｉｄｅｏ Ｄｅｎｏｉｓｉｎｇ ｂａｓｅｄ ｏｎ ａ ｓｐａｔｉｏｔｅｍｐｏｒａｌ Ｇａｕｓｓｉａｎ ｓｃａｌｅ ｍｉｘｔｕｒｅ ｍｏｄｅｌ」、ＩＥＥＥ Ｔｒａｎｓａｃｔｉｏｎｓ ｏｎ Ｃｉｒｃｕｉｔｓ ａｎｄ Ｓｙｓｔｅｍｓ ｆｏｒ Ｖｉｄｅｏ Ｔｅｃｈｎｏｌｏｇｙ、第２０巻、第７号、第１０３２−１０４０ページ、２０１０年において提示されるグローバルモーション補償（ＭＣ）方法を選択してもよい。それは簡潔で、高速で、信頼できる方法であり、整数画素の精度を提供する。Ｙ_３（ｍ，ｎ）は、すでにスケーリングおよび回転に関して補償されている補正された画像とし、Ｘ（ｍ，ｎ）を基準画像とする。また、２Ｄ ＦＦＴをＸおよびＹ_３に実施する。

そして、相互相関（ＣＣ）関数は以下のように定義される。

ここでＩＦＦＴ_２は逆２Ｄフーリエ変換、・は成分ごとの積、およびｃｏｎｊは複素共役を意味する。推定されるモーションベクトルは以下のように与えられる。

平行移動ベクトル推定モジュール１５２０において計算される、この推定されるモーションベクトルは、回転された歪画像Ｙ_３をレジストレーション画像Ｙ^＊としての最終位置に動かすために必要とされる、平行移動ベクトルと解釈される。すなわち、Ｙ^＊（ｍ，ｎ）＝Ｙ_３（ｍ−ｍ_ｏｐｔ，ｎ−ｎ_ｏｐｔ）。

図１８−図２２は、以下に記載する図２３において記載された本発明の第４の実施形態において提案されるレジストレーション方法の実行を表す実験の結果を示す。

図１８ａ−図１８ｅから図２２ａ−図２２ｅのそれぞれは、
−オリジナルの画像（Original）（ａ）、
−その歪んだバージョン（Distorted）（ｂ）、
−本発明の第３の実施形態の方法によってレジストレーションされた歪画像（Proposed）（ｃ）、
−ＳＵＲＦ方法によってレジストレーションされた歪画像（SURF）（ｄ）、および
−ＳＩＦＴ方法によってレジストレーションされた歪画像（SIFT）（ｅ）
を示す。

具体的に、図１８ａ、図１８ｂ、図１８ｃ、図１８ｄおよび図１８ｅはそれぞれ、オリジナルの画像、その歪んだバージョン（スケール＝０．８、回転＝０．２ラジアン、平行移動＝（４０，４０）、ノイズ標準偏差σ_ｎ＝１０）、本発明で提案の方法を用いたレジストレーション画像、ＳＵＲＦを用いたレジストレーション画像、ＳＩＦＴを用いたレジストレーション画像を示す。

図１９ａ、図１９ｂ、図１９ｃ、図１９ｄおよび図１９ｅはそれぞれ、オリジナルの画像、その歪んだバージョン（スケール＝０．８、回転＝０．２ラジアン、平行移動＝（４，４）、ノイズ標準偏差σ_ｎ＝５０）、本発明で提案の方法を用いたレジストレーション画像、ＳＵＲＦを用いたレジストレーション画像、ＳＩＦＴを用いたレジストレーション画像を示す。

図１８ａ−図１８ｅおよび図１９ａ−図１９ｅにおいて、３つの方法すべては、比較的小さいものからほどほどのノイズ（それぞれσ_ｎ＝１０およびσ_ｎ＝５０）と組み合わされた大きいものからほどほどの空間シフト（それぞれ４０，４０および４，４）による歪に対して良好に作用することがわかる。

図２０ａ、図２０ｂ、図２０ｃ、図２０ｄおよび図２０ｅはそれぞれ、オリジナルの画像、その歪んだバージョン（スケール＝０．８、回転＝０．２ラジアン、平行移動＝（４，４）、ノイズ標準偏差σ_ｎ＝１００）、本発明で提案の方法を用いたレジストレーション画像、ＳＵＲＦを用いたレジストレーション画像、ＳＩＦＴを用いたレジストレーション画像を示す。この場合には、ＳＵＲＦとＳＩＦＴはいずれも失敗であるが、提案の方法は良好に作用する。

図２１ａ、図２１ｂ、図２１ｃ、図２１ｄおよび図２１ｅはそれぞれ、オリジナルの画像、その歪んだバージョン（スケール＝０．５、回転＝０．２ラジアン、平行移動＝（４，４）、ノイズ標準偏差σ_ｎ＝１０）、本発明で提案の方法を用いたレジストレーション画像、ＳＵＲＦを用いたレジストレーション画像、ＳＩＦＴを用いたレジストレーション画像を示す。倍率が０．５では、３つの方法はいずれも良好に作用する。

図２２ａ、図２２ｂ、図２２ｃ、図２２ｄおよび図２２ｅはそれぞれ、オリジナルの画像、その歪んだバージョン（スケール＝０．１、回転＝０．２ラジアン、平行移動＝（４，４）、ノイズ標準偏差σ_ｎ＝１０）、提案の方法を用いたレジストレーション画像、ＳＵＲＦを用いたレジストレーション画像、ＳＩＦＴを用いたレジストレーション画像を示す。倍率が０．１のとき、本発明の実施形態で提案の方法は良好に作用する。しかしながら、ＳＵＲＦおよびＳＩＦＴはいずれもうまくいかない。このように、本発明の実施形態で提案の方法は、入力画像が非常に小さくスケーリングされた場合にＳＩＦＴおよびＳＵＲＦよりも良好である。

まとめると、σ_ｎ＝１００のとき、ＳＵＲＦおよびＳＩＦＴ方法はうまくいかないが、本発明で提案の方法はこの場合に良好に作用する。加えて、倍率ａ＝０．５のとき、３つの方法はすべて良好なレジストレーションを生成する（図２１を参照）。しかしながら、倍率ａ＝０．１のとき、ＳＵＲＦおよびＳＩＦＴはいずれもうまくいかないが、本発明で提案の方法は非常に良好に作用する（図２２ｃを参照）。

本発明の第３の実施形態において提案の方法の計算の複雑性が、２つの既存の方法と比較された。従来技術文献において、Ｂ．Ｓ．ＲｅｄｄｙおよびＢ．Ｎ．Ｃｈａｔｔｅｒｊｉ「Ａｎ ＦＦＴ−ｂａｓｅｄ ｔｅｃｈｎｉｑｕｅ ｆｏｒ ｔｒａｎｓｌａｔｉｏｎ，ｒｏｔａｔｉｏｎ ａｎｄ ｓｃａｌｅ ｉｎｖａｒｉａｎｔ ｉｍａｇｅ ｒｅｇｉｓｔａｔｉｏｎ」、ＩＥＥＥ Ｔｒａｎｓａｃｔｉｏｎｓ ｏｎ Ｉｍａｇｅ Ｐｒｏｃｅｓｓｉｎｇ、第５巻、第８号、第１２６６−１２７１ページ、１９９６年は、平行移動、回転、およびスケーリングのパラメータを推定する技術を提示する。彼らの方法は６つの順２Ｄ ＦＦＴおよび６つの逆２Ｄ ＦＦＴを使用し、時間を消費しエラーが生じやすい。本発明の第３の実施形態では、出願人は平均の比を用いて倍率を推定し、ラドン変換を用いて基準画像と歪画像の間の回転角を推定する。画像のラドン変換は、Ｙ．Ｐａｎ、Ｋ．ＬｉおよびＭ．Ｈａｍｄｉ「Ａｎ ｉｍｐｒｏｖｅｄ ｃｏｎｓｔａｎｔ−ｔｉｍｅ ａｌｇｏｒｉｔｈｍ ｆｏｒ ｃｏｍｐｕｔｉｎｇ ｔｈｅ ｒａｄｏｎ ａｎｄ ｈｏｕｇｈ ｔｒａｎｓｆｏｒｍｓ ｏｎ ａ ｒｅｃｏｎｆｉｇｕｒａｂｌｅ ｍｅｓｈ」、ＩＥＥＥ Ｔｒａｎｓａｃｔｉｏｎｓ ｏｎ Ｓｙｓｔｅｎｓ，Ｍａｎ， ａｎｄ Ｃｙｂｅｒｎｅｔｉｃｓ，Ｐａｒｔ Ａ−Ｓｙｓｔｅｍｓ ａｎｄ Ｈｕｍａｎｓ、第２９巻、第４号、第４１７−４２１ページ、１９９９年に示されるように、再構成可能なメッシュ上で一定の時間Ｏ（１）で実行可能である。Ｋ．Ｊａｆａｒｉ−ＫｈｏｕｚａｎｉおよびＨ．Ｓｏｌｔａｎｉａｎ−Ｚａｄｅｈ「Ｒａｄｏｎ ｔｒａｎｓｆｏｒｍ ｏｒｉｅｎｔａｔｉｏｎ ｅｓｔｉｍａｔｉｏｎ ｆｏｒ ｒｏｔａｔｉｏｎ ｉｎｖａｒｉａｎｔ ｔｅｘｔｕｒｅ ａｎａｌｙｓｉｓ」、ＩＥＥＥ Ｔｒａｎｓａｃｔｉｏｎｓ ｏｎ Ｐａｔｔｅｒｎ Ａｎａｌｙｓｉｓ ａｎｄ Ｍａｃｈｉｎｅ Ｉｎｔｅｌｌｉｇｅｎｃｅ、 第２７巻、第６号、第１００４−１００８ページ、２００５年に示されるように、多くとも、ラドン変換は２Ｄ ＦＦＴの複雑度で実施されることができる。そして出願人は２Ｄ ＦＦＴを利用してグローバルな空間シフトを得る。Ｇ．ＶａｒｇｈｅｓｅとＺ．Ｗａｎｇ「Ｖｉｄｅｏ ｄｅｎｏｉｓｉｎｇ ｂａｓｅｄ ｏｎ ａ ｓｐａｔｉｏｔｅｍｐｏｒａｌ Ｇａｕｓｓｉａｎ ｓｃａｌｅ ｍｉｘｔｕｒｅ ｍｏｄｅｌ」、ＩＥＥＥ Ｔｒａｎｓａｃｉｏｎｓ ｏｎ Ｃｉｒｃｕｉｔｓ ａｎｄ Ｓｙｓｔｅｍｓ ｆｏｒ Ｖｉｄｅｏ Ｔｅｃｈｎｏｌｏｇｙ、第２０巻、第７号、第１０３２−１０４０ページ、２０１０年を参照。本発明の第３の実施形態の方法は３つの順２Ｄ ＦＦＴ、２つの順２Ｄ ＦＦＴ、および２つのラドン変換のみを用いる。したがって、出願人の提示する方法は上に引用したＢ．Ｓ．ＲｅｄｄｙおよびＢ．Ｎ．Ｃｈａｔｔｅｒｊｉの論文に記載の従来技術の方法よりも高速である。Ｅ．Ｄｅ ＣａｓｔｒｏおよびＣ．Ｍｏｒａｎｄｉ「Ｒｅｇｉｓｔｒａｔｉｏｎ ｏｆ ｔｒａｎｓｌａｔｅｄ ａｎｄ ｒｏｔａｔｅｄ ｉｍａｇｅｄ ｕｓｉｎｇ ｆｉｎｉｔｅ Ｆｏｕｒｉｅｒ ｔｒａｎｓｆｏｒｍｓ」、ＩＥＥＥ Ｔｒａｎｓａｃｔｉｏｎｓ ｏｎ Ｐａｔｔｅｒｎ Ａｎａｌｙｓｉｓ ａｎｄ Ｍａｃｈｉｎｅ Ｉｎｔｅｌｌｉｇｅｎｃｅ、ｖｏｌ．ＰＡＭＩ−９５、第７００−７０３ページ、１９８７年において、著者は１８０の逆２Ｄ ＦＦＴを必要とする画像レジストレーション方法を提示する。そしてさらに、スケーリングの変化が存在する場合はその方法は失敗する。その逆に、本発明の第３の実施形態の方法はスケーリングの変化が存在する場合にも良好に作用し、Ｅ．Ｄｅ ＣａｓｔｒｏおよびＣ．Ｍｏｒａｎｄｉによって記載される方法よりも計算の複雑性がはるかに低い。

本発明の第３の実施形態の方法およびシステムは、既存の従来技術の方法に対して以下の利点を提供する。本発明の第３の実施形態の方法はノイズの多い環境でも非常に良好に作用するが、Ｄ．Ｇ．Ｌｏｗｅ「Ｄｉｓｔｉｎｃｔｉｖｅ Ｉｍａｇｅ Ｆｅａｔｕｒｅｓ ｆｒｏｍ Ｓｃａｌｅ−Ｉｎｖａｒｉａｎｔ Ｋｅｙｐｏｉｎｔｓ」、Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｊｏｕｒｎａｌ ｏｆ Ｃｏｍｐｕｔｅｒ Ｖｉｓｉｏｎ、第６０巻、第２号、第９１−１１０ページ、２００４年およびＨ．Ｂａｙ、Ａ．Ｅｓｓ、Ｔ．ＴｕｙｔｅｌａａｒｓおよびＬ．Ｖａｎ Ｇｏｏｌ「ＳＵＲＦ：Ｓｐｅｅｄｅｄ Ｕｐ Ｒｏｂｕｓｔ Ｆｅａｔｕｒｅｓ」、Ｃｏｍｐｕｔｅｒ Ｖｉｓｉｏｎ ａｎｄ Ｉｍａｇｅ Ｕｎｄｅｒｓｔａｎｄｉｎｇ （ＣＶＩＵ）、第１１０巻、第３号、第３４６−３５９ページ、２００８年によって提示されたものなどの既存の方法は、ノイズのレベルがあまりに高い場合はうまくいかない。加えて、本発明の実施形態の方法は、計算複雑性Ｏ（ＭＮｌｏｇ（ＭＮ）＋ＫＬｌｏｇ（Ｌ））に関する高速性をもつ。ここで、入力画像はＭ×Ｎの寸法であり、ラドン画像はＫ×Ｌの寸法である。さらに、本願発明の実施形態の方法は３つのパラメータのすべて（平行移動、回転、およびスケーリング）を推定することが可能である。これに対して大抵の既存の方法は、例えばＫ．Ｊａｆａｒｉ−ＫｈｏｕｚａｎｉおよびＨ．Ｓｏｌｔａｎｉａｎ−Ｚａｄｅｈ「Ｒａｄｏｎ ｔｒａｎｓｆｏｒｍ ｏｒｉｅｎｔａｔｉｏｎ ｅｓｔｉｍａｔｉｏｎ ｆｏｒ ｒｏｔａｔｉｏｎ ｉｎｖａｒｉａｎｔ ｔｅｘｔｕｒｅ ａｎａｌｙｓｉｓ」、ＩＥＥＥ Ｔｒａｎｓａｃｔｉｏｎｓ ｏｎ Ｐａｔｔｅｒｎ Ａｎａｌｙｓｉｓ ａｎｄ Ｍａｃｈｉｎｅ Ｉｎｔｅｌｌｉｇｅｎｃｅ、 第２７巻、第６号、第１００４−１００８ページ、２００５年、または、Ｗ．Ｗｅｉ、Ｓ．Ｗａｎｇ、Ｘ．ＺｈａｎｇおよびＺ．Ｔａｎｇ「Ｅｓｔｉｍａｔｉｏｎ ｏｆ ｉｍａｇｅ ｒｏｔａｔｉｏｎ ａｎｇｌｅ ｕｓｉｎｇ ｉｎｔｅｒｐｏｌａｔｉｏｎ−ｒｅｌａｔｅｄ ｓｐｅｃｔｒａｌ ｓｉｇｎａｔｕｒｅｓ ｗｉｔｈ ａｐｐｌｉｃａｔｉｏｎ ｔｏ ｂｌｉｎｄ ｄｅｔｅｃｔｉｏｎ ｏｆ ｉｍａｇｅ ｆｏｒｇｅｒｙ」、ＩＥＥＥ Ｔｒａｎｓａｃｔｉｏｎｓ ｏｎ Ｉｎｆｏｒｍａｔｉｏｎ Ｆｏｒｅｎｓｉｃｓ ａｎｄ Ｓｅｃｕｒｉｔｙ、第５巻、第３号、第５０７−５１７ページ、２０１０年において提示されたもののように、１つまたは２つのパラメータを計算するのみである。次の章において実施の実験結果は、特にノイズの多い環境において、出願人が提示した方法が画像のレジストレーションに適していることを示す。

本発明の第３の実施形態の実験結果
出願人は５１２×５１２のＢａｒｂａｒａ ａｎｄ Ｌｅｎａの画像においてある実験を行った。そして再びＬＩＶＥ画質評価データベースリリース２（Ｈ．Ｒ．Ｓｈｅｉｋｈ、Ｚ．Ｗａｎｇ、Ｌ．ＣｏｒｍａｃｋおよびＡ．Ｃ．Ｂｏｖｉｋ「ＬＩＶＥ ｉｍａｇｅ ｑｕａｌｉｔｙ ａｓｓｅｓｓｍｅｎｔ ｄａｔａｂａｓｅ ｒｅｌｅａｓｅ ２」、ｈｔｔｐ：／／ｌｉｖｅ．ｅｃｅ．ｕｔｅｘａｓ．ｅｄｕ／ｒｅｓｅａｒｃｈ／ｑｕａｌｉｔｙ）を使用した。このデータベースは、５つのタイプの歪みを用いて２９のオリジナルの画像から作成された７７９の歪画像を含む。歪みは、ＪＰＥＧ圧縮、ＪＰＥＧ２０００圧縮、ガウスホワイトノイズ（ＧＷＮ）、ガウスぼかし（ＧＢｌｕｒ）、およびレイリー高速フェーディング（ＦＦ）チャネルモデルを含む。出願人はＬＩＶＥ画像データベース内の７７９のすべての歪画像において実験を行った。出願人は本発明の第３の実施形態において提示した方法のためにＭＡＴＬＡＢコードを使用した。

表６および表７は、平行移動、回転、およびスケーリング歪みを有し、さらにノイズを加えた、Ｌｅｎａ ａｎｄ Ｂａｒｂａｒａの画像に関する実験的な結果を示す。いずれの表においても、出願人の推定するパラメータは、歪画像を生成した入力歪みパラメータに非常に近い。加えて、出願人の提示する方法はノイズの多い環境において非常に良好に作用する。これに対して、σｎ＝１００において、ＳＩＦＴはマッチングするキーポイントを発見することができず、ＳＵＲＦはＬｅｎａ画像に関して正確な結果を生成しない。

本発明の実施形態における提案の方法における回転角の精度が、点Ｌの数に依存することは、注目すべきである。出願人のシミュレーションにおいて、出願人はラドン変換の回転角においてＬ＝３６０のサンプル点を用いた。したがって精度は１度（０．０１７４ラジアンに等しい）である。サンプル点の数を増やせば正確性を上昇させることが可能だが、それにより計算の複雑性が増すことになる。したがって、パラメータの正確性と計算の複雑性の間にはトレードオフの関係がある。

表８はＬＩＶＥデータベース内のすべての歪画像が、倍率０．９、回転角０．１×１８０／π度、空間シフト水平方向に４画素、垂直方向に４画素で、さらに歪んだ場合の結果を表にしたものである。このデータベースに関して、ノイズが加えられない場合は、３つの方法（提案のもの、ＳＩＦＴ、およびＳＵＲＦ）すべてがパラメータの良好な推定を可能とする。出願人は、表６および表７に示すように、このデータベース内の画像にノイズを加えた場合、本発明の実施形態の提案の方法がＳＩＦＴおよびＳＵＲＦよりも良好に動作することを推定する。

表９は寸法５１２×５１２の１対の画像をレジストレーションする３つの方法に関する実行時間を秒で示した表である。本発明の実施形態で出願人が提示する方法はＳＵＲＦよりも消費時間が短いが、ＳＩＦＴよりは長いことがわかる。ＳＩＦＴはＣ言語で記述されており、したがってＭＡＴＬＡＢよりも実行が速いことに留意する。Ｈ．Ｂａｙ、Ａ．Ｅｓｓ、Ｔ．ＴｕｙｔｅｌａａｒｓおよびＬ．Ｖａｎ Ｇｏｏｌ「ＳＵＲＦ：Ｓｐｅｅｄｅｄ Ｕｐ Ｒｏｂｕｓｔ Ｆｅａｔｕｒｅｓ」、Ｃｏｍｐｕｔｅｒ Ｖｉｓｉｏｎ ａｎｄ Ｉｍａｇｅ Ｕｎｄｅｒｓｔａｎｄｉｎｇ （ＣＶＩＵ）、第１１０巻、第３号、第３４６−３５９ページ、２００８年によれば、仮にＳＩＦＴおよびＳＵＲＦが共にＭＡＴＬＡＢで記述されていたならば、ＳＩＦＴはＳＵＲＦほど速くない。したがって、本発明の実施形態で出願人が提示する方法はＳＵＲＦおよびＳＩＦＴのいずれよりも速いと結論づけることができる。加えて、出願人が実験において示したように、出願人が提示する方法はノイズが多い環境においてＳＵＲＦおよびＳＩＦＴよりも良好に動作する。

表６。平行移動された（４，４）、回転された（０．１ラジアン）、およびスケーリングされた（０．９）Ｌｅｎａ画像を、異なるノイズレベルでマッチングした結果である。

表７。平行移動された（４，４）、回転された（０．１ラジアン）、およびスケーリングされた（０．９）Ｂａｒｂａｒａ画像を、異なるノイズレベルでマッチングした結果である。

表８。ＬＩＶＥ画像データベース内のすべてのすでに歪画像に関して、倍率０．９、回転角０．１×１８０／π度、空間シフト水平方向に４画素、垂直方向に４画素でさらに歪ませたものの、全体として推定されるパラメータを示す。

表９。サイズ５１２×５１２の１対の画像をレジストレーションのための実行時間を示す。ＳＵＲＦおよび提案の方法はＭＡＴＬＡＢコードを用いるがＳＩＦＴはＣコードを用いる。Ｈ．Ｂａｙ、Ａ．Ｅｓｓ、Ｔ．ＴｕｙｔｅｌａａｒｓおよびＬ．Ｖａｎ Ｇｏｏｌ「ＳＵＲＦ：Ｓｐｅｅｄｅｄ Ｕｐ Ｒｏｂｕｓｔ Ｆｅａｔｕｒｅｓ」、Ｃｏｍｐｕｔｅｒ Ｖｉｓｉｏｎ ａｎｄ Ｉｍａｇｅ Ｕｎｄｅｒｓｔａｎｄｉｎｇ （ＣＶＩＵ）、第１１０巻、第３号、第３４６−３５９ページ、２００８年によれば、仮にＳＩＦＴおよびＳＵＲＦが共にＭＡＴＬＡＢを用いたならば、ＳＩＦＴはＳＵＲＦよりも遅い。

本発明の第３の実施形態の代替例が、以下に記載される第４の実施形態において開示される。

画像前処理システム１３０の第４の実施形態は第３の実施形態の拡張バージョンであり、歪画像の回転およびサイズ変更の順番が異なっていることを特徴とする。第３の実施形態の動作を要約する。歪画像Ｙがまずセンタリングされてセンタリングされた画像Ｙ_１となる。そしてサイズ変更されて（スケーリングされて）、基準画像とサイズをマッチングされて、正規化された画像Ｙ_２となる。回転角は基準画像Ｘおよび正規化された画像Ｙ_２のラドン変換の相関角によって決定される。そして正規化された画像Ｙ_２は回転されて、基準画像Ｘと方向をアライメントされ、回転された歪画像Ｙ_３となる。そして最後に、回転された歪画像Ｙ_３は水平方向および垂直方向に平行移動されて基準画像Ｘとアライメントされ、それによってレジストレーション画像Ｙ^＊となる。
第４の実施形態は、ラドン変換を計算する前に歪画像Ｙがスケーリングされるのではなく、第４の実施形態においてはラドン変換を計算する前に基準画像Ｘが歪画像Ｙとサイズを同じにされて、歪画像の必要なサイズ変更は回転および平行移動が行われた後で、最後に行われる点で、顕著に異なる。

記載の読解を容易にするために、第３の実施形態のいくつかの変数名は第４の実施形態においても再び使用することに留意する。

図２３は、図１の画像前処理システム１３０の第４の実施形態である拡張ラドン変換に基づく画像アライメントシステム２３００のブロック図であり、以下のモジュールを含む。
倍率「ａ」を決定する、図１５の倍率推定モジュール１５０２。

センタリングされた画像Ｙ_１を生成する、オプションのセンタリングモジュール２３０５。それは図１５の画像マスキングモジュール１５０４、重心計算モジュール１５０６、および画像センタリングモジュール１５０８を含む。

倍率「ａ」の逆数によって基準画像Ｘをスケーリングして、サイズ変更された基準画像Ｘ_２とする、画像前スケーリングモジュール２３１０。

サイズ変更された基準画像Ｘ_２およびセンタリングされた画像Ｙ_１それぞれのラドン変換Ｒ_１およびＲ_２を生成するラドン変換モジュール１５１２、および、ラドン変換Ｒ１およびＲ２から角度を抽出して、角度からシフト効果を消滅させて、回転角「θ_０」を決定する、回転角推定ブロック２３２０を含む第２の回転角決定ユニット２３１５。

センタリングされた画像Ｙ_１が回転角「θ_０」によって回転されて、回転された歪画像Ｙ_３を形成するが、回転された歪画像Ｙ_３はまだ基準画像Ｘから横方向にオフセットされている場合のある、図１５の画像回転モジュール１５１８。

回転された歪画像Ｙ_３とサイズ変更された基準画像Ｘ_２の間のオフセットベクトル「ＴＶ」が決定される、図１５の平行移動推定モジュール１５２０。

オフセットベクトル「ＴＶ」が、回転された歪画像Ｙ_３に適用され、補償された画像Ｙ_４が生成される、図１５の画像平行移動モジュール１５２２。

補償された画像Ｙ_４を倍率「ａ」で基準画像Ｘのサイズにスケーリングし、そして、レジストレーション画像Ｙ^＊を生成する、画像スケーリングモジュール２３３０。

拡張ラドン変換に基づく画像アライメントシステム２３００に用いられるモジュール１５０２、１５１２、１５１８、１５２０および１５２２は、ここでは、図１５の第１のラドン変換に基づく画像アライメントシステム１５００内のものと同じ機能を有するが、それらの入力データおよび出力データは変更されることがある。

拡張ラドン変換に基づく画像アライメントシステム２３００の、画像Ｘの前スケーリングの役割を果たすモジュール２３１０と、画像Ｙ４の最後のスケーリングを行ってレジストレーション画像Ｙ^＊を生成する２３３０は、変更されたトポロジーに適応している。

図１５の画像マスキングモジュール１５０４、重心計算モジュール１５０６、および画像センタリングモジュール１５０８を含み、同様の機能を有するモジュール２３０５は、拡張ラドン変換に基づく画像アライメントシステム２３００ではオプションである。画像センタリングモジュールの機能は、画像Ｘを前スケーリングする前に画像Ｘをセンタリングするために用いられてもよい。

拡張ラドン変換に基づく画像アライメントシステム２３００の変形において、図１４の構成ブロックは異なって組み合わされてもよい。例えば、回転角を決定する目的で、前スケーリングされた画像Ｘの代わりにセンタリングされた画像Ｙ_１を前スケーリングすることも可能である。

図２４は、拡張ラドン変換に基づく画像アライメントシステム２３００のモジュールにおいて例示することが可能な、第２のラドン変換に基づく方法２４００を示す。それは以下のステップを含む。

２４１０「倍率を推定する」
２４２０「マスクを生成する」
２４３０「画像をセンタリングする」
２４４０「基準画像Ｘをサイズ変更して、サイズ変更された基準画像Ｘ_２を生成する」
２４５０「回転角を推定する」
２４６０「歪画像を回転する」
２４７０「平行移動を推定する」
２４８０「画像を平行移動する」および
２４９０「歪画像をサイズ変更する」。

基準画像Ｘおよび歪画像Ｙが同じ数の画素を有していない場合、前記ステップ３１０「ゼロ埋込画像」において小さい方の画像の周囲にゼロを埋め込む。ここからは、ＸおよびＹはゼロが埋め込まれたバージョンのＸおよびＹであるとする。

図１６のステップ１６１０と同等のステップ２４１０「倍率を推定する」において、前記したような倍率「ａ」が決定される。

歪画像Ｙをセンタリングするためのステップ２４２０「マスクを生成する」、および２４３０「画像をセンタリングする」（それぞれ、第１のラドン変換に基づく方法１６００の図１６のステップ１６２０および１６３０と同等）は、第２のラドン変換に基づく方法２４００においてはオプションであり、ここでは必要とされない場合がある。なぜならより堅牢なラドン変換操作が用いられるからである。しかしながら、画像の実際の回転が多くの画素を失わせる（フレームの外側に位置させる）結果となる場合、パラメータの推定は不正確であり得る。したがって、堅牢性を高めるため、重心を計算し、画像をアライメントして（ステップ１６２０および１６３０）、歪画像Ｙをセンタリングされた画像Ｙ_１に変換することを選択してもよい。同様に、ここには示さないが、基準画像Ｘをオプションとしてセンタリングしてもよい。

ステップ２４４０「基準画像Ｘをサイズ変更して、サイズ変更された基準画像Ｘ_２を生成する」において、計算された倍率「ａ」の逆数１／ａが用いられ、基準画像が歪画像Ｙと同じスケールとなるように、基準画像Ｘからサイズ変更された基準画像Ｘ_２を生成する。このステップは回転角の推定において有利である。サイズ変更された基準画像をＸ_２（ｍ，ｎ）＝Ｘ（ｍ＊ａ，ｎ＊ａ）と表すものとする。

ステップ２４５０「回転角を推定する」において、サイズ変更された基準画像Ｘ_２および歪画像Ｙ_１のそれぞれのラドン変換Ｒ１およびＲ２の出力の相関から、基準画像Ｘおよび歪画像Ｙ_１の間の回転角「θ_０」を推定する。

図２５はステップ２４５０「回転角を推定する」の展開であり、以下のステップを有する。
２５１０「堅牢なラドン変換を計算する」
２５２０「Ｒ１およびＲ２の１Ｄ順ＦＦＴを計算してＲ３を生成する」
２５３０「１Ｄ逆ＦＦＴをＲ３に適用してＲ４を生成する」
２５４０「Ｒ４の各行において最大値をとる指数（インデックス）を計算する」および
２５５０「最大値の中央値からθ_０を計算する」。

ステップ２５１０「堅牢なラドン変換Ｒ１およびＲ２を計算する」において、ラドン変換Ｒ１およびＲ２はそれぞれ、サイズ変更された基準画像Ｘ_２およびセンタリングされた画像Ｙ_１（ステップ１６２０および１６３０における画像のセンタリングがスキップされた場合には、代わりに歪画像Ｙ）に関して計算される。２Ｄ離散関数Ａ（ｘ，ｙ）のラドン変換は、Ａ．Ａｖｅｒｂｕｃｈ、Ｒ．Ｒ．Ｃｏｉｆｍａｎ、Ｄ．Ｌ．Ｄｏｎｏｈｏ、Ｍ．Ｉｓｒａｅｌｉ、Ｙ．Ｓｈｋｏｌｎｉｓｋｙ、およびＩ．ＳｅｄｅｌｎｉｋｏｖのＡ Ｆｒａｍｅｗｏｒｋ ｆｏｒ Ｄｉｓｃｒｅｔｅ Ｉｎｔｅｇｒａｌ Ｔｒａｎｓｆｏｒｍａｔｉｏｎｓ ＩＩ−Ｔｈｅ ２Ｄ Ｄｉｓｃｒｅｔｅ Ｒａｄｏｎ Ｔｒａｎｓｆｏｒｍ、ＳＩＡＭ Ｊｏｕｒｎａｌ ｏｎ Ｓｃｉｅｎｔｉｆｉｃ Ｃｏｍｐｕｔｉｎｇ、３０（２）、第７８５−８０３ページ、２００８年に従い、方程式（２９）で定義される。

ここでδ（）はＤｉｒａｃのデルタ関数である。

Ａｖｅｒｂｕｃｈ ｅｔ ａｌ．によれば、Ｘ_２およびＹ（センタリングがスキップされなかった場合はＹ_１）の離散ラドン変換を実行可能であり、それぞれＲ１（ｒ，θ）、Ｒ２（ｒ，θ）と表され、両者とも寸法はＫ×Ｌである。

そしてＲ１（ｒ，θ）とＲ２（ｒ，θ）の両者にｒ方向に沿って１Ｄ ＦＦＴを行い、得られる行列の絶対値を得ることができる。これは角θの直線に沿うシフトを消去する。例えば、回転の位置が（ｘ０、ｙ０）の場合、角θの直線の空間的シフトは
ｂ＝ｘ０ｃｏｓ（θ）＋ｙ０ｃｏｓ（θ）
である。

この直線に沿ってＦＦＴを行い、大きさを得ることにより、後にシフトｂを消去することができる。このことは回転角とグローバルな空間シフトの両者の正確性を向上させる。

２つの画像（サイズ変更された基準画像Ｘ_２および歪画像ＹまたはそれのセンタリングされたバージョンＹ_１）が同一であるか、それらの間の回転およびありえるノイズを除いて少なくとも非常に類似しているという前提に基づくと、それらのラドン変換Ｒ_１（ｒ，θ）とＲ_２（ｒ，θ）は、θ方向に沿う回転シフトを除いて実質的に同一であることがわかる。すなわちＲ_２（ｒ，θ）＝Ｒ_１（ｒ，θ＋θ_０）である。θ_０を計算するために回転相互相関を用いることができるが、これはＯ（ＫＬ^２）の計算的複雑性で時間を消費する。なぜなら、Ｒ_１とＲ_２のすべての行に関して、Ｒ_２を回転シフトして、２つの行の間の相互相関を計算する必要があるからである。これはＯ（Ｌ^２）の計算的複雑性がある。全部でＫ個の行があるから、全部の計算的複雑性はＯ（ＫＬ^２）となるであろう。

より効率的な代替例として、代わりに１次元（１Ｄ）ＦＦＴを用いることができ、それはＯ（ＫＬｌｏｇＬ）の計算的複雑性がある。

第４の実施形態の好ましい方法が基礎を置くＦＦＴによる相互相関を計算する方法は、第３の実施形態において使用されたものと同一であり、上記ステップ１６６０の全体的記載を参照すべきであるが、上記シフト効果の消去により増大する。

上記考察に基づき、高速相互相関を用いてθ_０を求める。それはＯ（ＫＬｌｏｇＬ）の計算的複雑性がある。

Ｒ_１およびＲ_２のｍ番目の行をそれぞれｒ_１（ｎ）およびｒ_２（ｎ）と表すものとする。ここでｎ∈［１，Ｌ］である。相互相関が最大値になるオフセットが探される。すなわち、２つの行ｒ_１（ｎ）およびｒ_２（ｎ）が最もマッチングするオフセットが探される。

ステップ２５２０「Ｒ１およびＲ２の１Ｄ順ＦＦＴを計算してＲ３を生成する」において、１Ｄ順ＦＦＴがｒ_１（ｎ）およびｒ_２（ｎ）のそれぞれで実行され、第１の中間結果Ｒ３＝ｒ_３（ｕ）を生成する。

ここで「ｃｏｎｊ」は複素共役を表し、・は成分ごとの積である。

ステップ２５３０「１Ｄ逆ＦＦＴをＲ３に適用してＲ４を生成する」において、１Ｄ逆ＦＦＴがｒ_３（ｕ）に適用される。すなわち、第２の中間結果Ｒ４＝ｒ_４（ｎ）＝ＩＦＦＴ（ｒ_３（ｕ））が計算される。第２の中間結果Ｒ４は、ラドン変換Ｒ１およびＲ２の行の組み合わせの間の相互相関値のベクトルを構成する。
ステップ２５４０「Ｒ４の各行の中で最大値をとる指数を計算する」において、すべての行に対して、得られた第２の中間結果Ｒ４＝ｒ_４（ｎ）の最大値の位置指数θ^＊（ｍ）を求める。すなわち

である。

ステップ２５５０「最大値の中央値からθ_０を計算する」において、θ_０＝ｍｅｄｉａｎ（θ^＊（ｍ））として、これらの最大値の中から中央値を取る。
以下のＭＡＴＬＡＢコードサンプル３が回転角の推定を実施する。

ここで図２４の記載に戻る。

ステップ２４６０「歪画像を回転する」において、歪画像Ｙは方向の違いを補償するように−θ度だけ回転される。回転の結果が回転された歪画像Ｙ_３である。

ステップ２４７０「平行移動を推定する」において、サイズ変更された基準画像Ｘ_２にアライメントするために必要な、回転された歪画像Ｙ_３の平行移動が、Ｘ_２およびＸ_３の関数として決定される。

平行移動のために、Ｇ．ＶａｒｇｈｅｓｅとＺ．Ｗａｎｇ「Ｖｉｄｅｏ ｄｅｎｏｉｓｉｎｇ ｂａｓｅｄ ｏｎ ａ ｓｐａｔｉｏｔｅｍｐｏｒａｌ Ｇａｕｓｓｉａｎ ｓｃａｌｅ ｍｉｘｔｕｒｅ ｍｏｄｅｌ」、ＩＥＥＥ Ｔｒａｎｓａｃｉｏｎｓ ｏｎ Ｃｉｒｃｕｉｔｓ ａｎｄ Ｓｙｓｔｅｍｓ ｆｏｒ Ｖｉｄｅｏ Ｔｅｃｈｎｏｌｏｇｙ、第２０巻、第７号、第１０３２−１０４０ページ、２０１０年に記載のグローバルモーション補償（ＭＣ）方法を選択してもよい。それは簡潔で、高速で、信頼できる方法であり、整数画像の精度を提供する。Ｙ_３（ｍ，ｎ）はすでに回転の補償をされた画像であるとし、Ｘ_２（ｍ，ｍ）はスケーリングされた基準画像であるとする。また、

である。
すると、高速相互相関（ＦＣ）関数を方程式（３２）において定義できる。

ここで、ＩＦＦＴ_２は逆２Ｄフーリエ変換であり、ｃｏｎｊは複素共役であり、記号・は成分ごとの積を示す。推定されたモーションベクトル（平行移動ベクトル）は方程式（３３）で与えられる。

ステップ２４８０「画像を平行移動する」において、計算された平行移動ベクトルは回転された歪画像Ｙ_３に適用されて平行移動を補償し、平行移動された画像Ｙ_４をもたらす。この関数を形式的にＹ_４（ｍ，ｎ）＝Ｙ_３（ｍ−ｍ_ｏｐｔ，ｎ−ｎ_ｏｐｔ）と示す。

ステップ２４９０「歪画像をサイズ変更する」において、最終的な画像、すなわちレジストレーション画像Ｙ^＊を得るために、歪画像Ｙ４が基準画像Ｘと同じスケールを有するように、計算された倍率ａが用いられて歪画像Ｙ４がサイズ変更される。それは今や平行移動、回転、およびスケーリングに関して補償されたものである。この関数を形式的にＹ^＊（ｍ，ｎ）＝Ｙ_４（ｍ／ａ，ｎ／ａ）と示す。図２３に示す本発明の第４の実施形態のシステムは、部分的にラドン変換のために用いられる処理であるため非常に堅牢であり、第３の実施形態のシステムとは対照的に、いかなる回転角に対しても良好に作用する。

図２６は、確認の目的で実施される画像アライメントシステム１００の実施形態の組み合わせられたブロック図２６００を示す。それは、ＣＰＵ２６０４、ネットワークＩ／Ｏシステム、およびコマンドインターフェイス２６０８を含むプロセッサ２６０２、メモリ、ＤＶＤ、ＣＤ−ＲＯＭなどの不揮発性コンピュータ読み取り可能記憶媒体の形態の、前記プロセッサによる実行のためのコンピュータ読み取り可能な命令が記憶されたコンピュータメモリ２６１０、および、前記画像アライメントシステム１００による処理のための基準画像および試験画像を含む画像データベース２６１２、を含む。前記コンピュータメモリ２６１０に記憶された前記コンピュータ読み取り可能な命令は、本発明の実施形態の少なくとも１つを実施するためのソフトウェアモジュールを含む。具体的には、アフィン変換に基づく画像アライメントシステム２００．１、代替アフィン変換に基づく画像アライメントシステム２００．２、第１のラドン変換に基づく画像アライメントシステム１５００、拡張ラドン変換に基づく画像アライメントシステム２３００である。コンピュータメモリ２６１０は本発明の実施形態における実施の際に用いられる、画像評価処理１６０などの画像評価プログラム２６１４、および、例えばＭＡＴＬＡＢ環境２６１８などのソフトウェアライブラリ２６１４のストレージをさらに有してもよい。

図２６は、本発明の異なる実施形態の個々のモジュールを、重複するモジュールの組として示している。例えばアフィン変換に基づく画像アライメントシステム２００．１と代替アフィン変換に基づく画像アライメントシステム２００．２はモジュール２３０、２４０、２５０、２６０および２７０を共用する。共用するモジュールに加えて、アフィン変換に基づく画像アライメントシステム２００．１はモジュール２２０および２８０をさらに含み、代替アフィン変換に基づく画像アライメントシステム２００．２はモジュール２９０を含む。同様に、第１のラドン変換に基づく画像アライメントシステム１５００および拡張ラドン変換に基づく画像アライメントシステム２３００はモジュールの組を共用する。すなわち１５０２、１５１２、１５１８、１５２０および１５２２である。共用するモジュールに加えて、第１のラドン変換に基づく画像アライメントシステム１５００はさらにモジュール１５０４、１５０６、１５０８、１５１０および１５１６を用い、拡張ラドン変換に基づく画像アライメントシステム２３００は２３０５、２３１０、２３２０および２３３０を含む。

本発明の実施形態は、マルチコアＣＰＵであってもよいＣＰＵと、コンピュータ読み取り可能な媒体、例えばメモリ、ＤＶＤ、ＣＤ−ＲＯＭ、フロッピー（登録商標）、磁気テープ、またはその他の記憶媒体であって、前記ＣＰＵによって実行された場合に、上記のようなシステムのモジュールを形成するコンピュータ読み取り可能な命令が記憶されている媒体と、を有する、汎用または専用コンピュータを含む。代替として、図２ａ、図２ｂ、図１５、および図２３は、上記したこれらのシステムのモジュールを形成するために前記ＣＰＵによって実行されるコンピュータ読み取り可能な命令を有するコンピュータ読み取り可能な記憶媒体を有する専用コンピュータとファームウェアとの組み合わせまたは特別な専用のハードウェアを含んでもよい。図２ａ、図２ｂ、図１５、および図２３のシステムの各モジュールは、ファームウェアまたは、代替として、図２６に示すようにプロセッサによって実行するための、コンピュータ読み取り可能な記憶媒体内に記憶されたコンピュータ読み取り可能な命令を含んでもよい。

結論として、本発明は歪画像を前処理して基準画像にアライメントするための改良されたシステムおよび方法を提供し、それらはより正確な視覚的な画質評価を提供するように使用されてもよく、同様に、歪画像を基準画像にレジストレーションするように作用するように使用されてもよい。

本発明の実施形態は詳細に記載されているが、実施形態の変形や変更を以下の特許請求の範囲内で行い得ることは当業者には明らかである。
下記は、本願に当初より記載の発明である。
＜請求項１＞
歪画像Ｙを基準画像Ｘにアライメントされたレジストレーション画像Ｙ^＊に処理する画像レジストレーション方法であって、
（ａ）前記基準画像Ｘと前記歪画像Ｙとの間の倍率「ａ」を決定するステップと、
（ｂ）前記基準画像Ｘを前記倍率「ａ」の逆数でサイズ変更し、それによって、サイズ変更された基準画像Ｘ２を生成するステップと、
（ｃ）前記サイズ変更された基準画像Ｘ２と前記歪画像Ｙの間の回転角「θ_０」を決定するステップと、
（ｄ）前記歪画像Ｙを回転角「−θ_０」だけ回転することによって、回転された歪画像Ｙ３を決定するステップと、
（ｅ）前記回転された歪画像Ｙ３を前記倍率「ａ」でサイズ変更し、それによって、前記レジストレーション画像Ｙ^＊を得るステップと、を含む、
画像レジストレーション方法。
＜請求項２＞
前記ステップ（ａ）は、前記基準画像Ｘの画素値の和と前記歪画像Ｙの画素値の和の比を計算することにより倍率「ａ」を決定するステップを含む、請求項１に記載の方法。
＜請求項３＞
前記ステップ（ａ）は、前記基準画像Ｘと前記歪画像Ｙのうち小さい方をゼロ値の画素で埋込み、前記基準画像Ｘと前記歪画像Ｙの水平方向および垂直方向の寸法（ｍおよびｎ）を等しくするステップをさらに含む、請求項１または２に記載の方法。
＜請求項４＞
前記ステップ（ｃ）は、前記歪画像Ｙをセンタリングすることでセンタリングされた歪画像Ｙ１を生成するステップと、前記サイズ変更された基準画像Ｘ２と前記センタリングされた歪画像Ｙ１の間の回転角「θ_０」の決定を行うステップと、をさらに含む、請求項１から３のいずれか一項に記載の方法。
＜請求項５＞
前記ステップ（ｃ）は、前記サイズ変更された基準画像Ｘ２および前記歪画像Ｙのラドン変換Ｒ１およびＲ２をそれぞれ形成するステップと、記回転角「θ_０」を決定するために前記ラドン変換Ｒ１およびＲ２を用いるステップと、を含む、請求項１から４のいずれか一項に記載の方法。
＜請求項６＞
前記ステップ（ｃ）は、前記サイズ変更された基準画像Ｘ２と前記センタリングされた歪画像Ｙ１のラドン変換Ｒ１およびＲ２をそれぞれ形成するステップと、前記ラドン変換Ｒ１およびＲ２の行の相互相関を計算して回転角「θ_０」を決定するステップと、を含む、請求項４または５に記載の方法。
＜請求項７＞
前記ラドン変換Ｒ１およびＲ２の行の相互相関を計算して回転角「θ_０」を決定するステップは、
（ｉ）前記ラドン変換Ｒ１の各行と前記ラドン変換Ｒ２の各行の間の円周方向相互相関の組を計算するステップであって、前記各円周方向相互相関は前記各行の間の回転オフセット角「θ」を規定するものである、該ステップと、
（ｉｉ）各行に対して、最も高い値を有する円周方向相互相関を選択するステップと、
（ｉｉｉ）各行に対して各選択された円周方向相互相関によって規定される回転オフセット「θ」を決定し、回転角「θ_０」を、前記決定された回転オフセット「θ」の中央値に等しくなるよう設定するステップと、をさらに含む、
請求項５または６に記載の方法。
＜請求項８＞
前記レジストレーション画像Ｙ^＊の視覚品質評価をさらに実行する、請求項１から７のいずれか一項に記載の方法。
＜請求項９＞
歪画像Ｙを基準画像Ｘにアライメントすることによってレジストレーション画像Ｙ^＊に処理するシステムであって、
プロセッサと、
前記プロセッサによって実行するためのコンピュータ読み取り可能な命令を有するメモリ装置と、を含み、
前記基準画像Ｘと前記歪画像Ｙの間の倍率「ａ」を決定する倍率推定モジュールと、
前記基準画像Ｘと前記歪画像Ｙの間の回転角「θ_０」を推定する回転角決定モジュールと、
前記歪画像Ｙを回転角「−θ_０」だけ回転することによって回転された歪画像Ｙ３を形成する画像回転モジュールと、
前記回転された歪画像Ｙ３をサイズ変更してレジストレーション画像Ｙ^＊を生成する画像スケーリングモジュールと、
を形成するシステム。
＜請求項１０＞
センタリングされた画像Ｙ１を生成するための任意的なセンタリングモジュールをさらに含み、前記回転角決定モジュールは、前記基準画像Ｘと前記センタリングされた画像Ｙ１の間の回転角「θ_０」を推定するように構成される、請求項９に記載のシステム。
＜請求項１１＞
前記倍率「ａ」の逆数で前記基準画像Ｘをスケーリングしてサイズ変更された基準画像Ｘ２とする画像前スケーリングモジュールをさらに含み、前記回転角決定モジュールは、前記サイズ変更された基準画像Ｘ２と前記センタリングされた画像Ｙ１の間の回転角「θ_０」を推定するように構成される、請求項９または１０に記載のシステム。
＜請求項１２＞
前記回転角決定モジュールはさらに、前記サイズ変更された基準画像Ｘ２と前記歪画像Ｙのラドン変換Ｒ１およびＲ２をそれぞれ形成し、回転角「θ_０」を決定するために前記ラドン変換Ｒ１およびＲ２を用いるように構成される、請求項１１に記載のシステム。
＜請求項１３＞
前記回転された歪画像Ｙ３と前記基準画像Ｘの間のオフセットベクトル「ＴＶ」を決定する平行移動推定モジュールと、
前記回転された歪画像Ｙ３を前記オフセットベクトル「ＴＶ」を用いて平行移動させ、レジストレーション画像Ｙ^＊を生成するように構成される画像平行移動モジュールと、
をさらに含む、請求項９から１２のいずれか一項に記載のシステム。
＜請求項１４＞
レジストレーション画像Ｙ^＊の視覚品質評価を実行するための画像評価モジュール処理モジュールをさらに有する、請求項９から１３のいずれか一項に記載のシステム。
＜請求項１５＞
前記画像評価モジュール処理モジュールは、前記視覚品質評価を、ピーク信号ノイズ比（ＰＳＮＲ）の決定、構造的類似性（ＳＳＩＭ）指標の計算および視覚情報忠実度（ＶＩＦ）指標の計算のいずれかによって実行するように構成される、請求項１４に記載のシステム。