(57)【要約】  本発明は、画像レジストレーション方法に関し、具体的には、画像または任意の次元の信号の間の並進シフトを識別する画像レジストレーション方法に関する。本明細書で説明する複数の画像をレジスタする方法は、１つまたは複数の空間領域または周波数領域の関数、フィルタ、またはカーネルの周波数領域表現を入手することであって、各関数は、少なくとも１つの画像特性を強調する、入手することと、関数を単一の周波数領域表現関数に組み合わせることと、相対変位を判定するために、組み合わされた周波数領域表現関数を画像に適用することとを含む。別の態様では、複数の画像をレジスタする方法は、画像の間のシフトを計算するために関数を適用することと、位相データ異常値を除外するために位相アンラッピング技法を適用することと、画像の間の副画素シフトを計算するために位相データを使用することとをも含む。【選択図】図１
2019517073.tif 000002
【特許請求の範囲】
【請求項１】
複数の画像のレジストレーションの方法であって、
１つまたは複数の空間領域または周波数領域の関数、フィルタ、またはカーネルの周波数領域表現を入手することであって、各関数は、少なくとも１つの画像特性を強調する、入手することと、
関数を単一の周波数領域表現関数に組み合わせることと、
相対変位を判定するために、前記組み合わされた周波数領域表現関数を前記画像に適用することと
を含む方法。
【請求項２】
前記複数の画像は、同一のおよび／または異なる次元および／またはサイズを有する、請求項１に記載の複数の画像のレジストレーションの方法。
【請求項３】
前記関数は、前記空間領域で定義される平滑化関数とすることができる、請求項１または２に記載の複数の画像のレジストレーションの方法。
【請求項４】
前記関数によって強調される前記画像特性のうちの１つは、雑音の低減された影響および／または画像信号対雑音比の向上である、請求項１から３のいずれか１項に記載の複数の画像のレジストレーションの方法。
【請求項５】
前記関数は、前記画像への適用の前に事前に計算され得る、請求項１から４のいずれか１項に記載の複数の画像のレジストレーションの方法。
【請求項６】
前記方法は、相互相関の周波数領域表現に前記関数の前記周波数領域表現を適用することをさらに含む、請求項１から５のいずれか１項に記載の複数の画像のレジストレーションの方法。
【請求項７】
前記方法は、事前に計算された周波数領域関数のライブラリをさらに含む、請求項１から６のいずれか１項に記載の複数の画像のレジストレーションの方法。
【請求項８】
前記周波数領域関数は、画像データ・セットに依存して前記ライブラリから選択される、請求項１から７のいずれか１項に記載の複数の画像のレジストレーションの方法。
【請求項９】
画像の間のシフトを計算するために関数を適用することと、
位相データ異常値を除外するために位相アンラッピング技法を適用することと、
前記画像の間の副画素シフトを計算するために位相データを使用することと
を含む、複数の画像をレジスタする方法。
【請求項１０】
整数部分と副画素部分との両方でのシフトの前記計算は、同一のおよび／または異なる次元および／またはサイズを有する前記複数の画像の間で行われる、請求項９に記載の複数の画像をレジスタする方法。
【請求項１１】
前記画像の間のシフトの前記計算は、１つまたは複数の関数の周波数領域表現を入手するステップをさらに含む、請求項９から１０のいずれか１項に記載の複数の画像をレジスタする方法。
【請求項１２】
前記関数は、前記周波数領域で定義される、請求項９から１１のいずれか１項に記載の複数の画像をレジスタする方法。
【請求項１３】
前記副画素部分でのシフトの前記計算は、前記位相データの勾配を推定するステップをさらに含む、請求項９から１２のいずれか１項に記載の複数の画像をレジスタする方法。
【請求項１４】
位相勾配の前記推定は、位相誤差の推定値および／または位相成分をアンラッピングすることおよび／または位相異常値を除去することに基づいて位相雑音を低減する方法を含む、請求項９から１３のいずれか１項に記載の複数の画像をレジスタする方法。
【請求項１５】
位相アンラッピング成分は、不適切にラップされた位相データを識別する、請求項９から１４のいずれか１項に記載の複数の画像をレジスタする方法。
【請求項１６】
誤差アンラッピング成分は、より正確な勾配推定値を入手するために、不適切にラップされた位相データをアンラップするための位相誤差の推定値を含む、請求項９から１５のいずれか１項に記載の複数の画像をレジスタする方法。
【請求項１７】
前記関数は、雑音低減プロパティを有する、請求項９から１６のいずれか１項に記載の複数の画像をレジスタする方法。
【請求項１８】
前記関数は、平滑化微分器関数をさらに含むことができる、請求項９から１７のいずれか１項に記載の複数の画像をレジスタする方法。
【請求項１９】
前記関数は、少なくとも１つの形状保存特性を含む、請求項９から１８のいずれか１項に記載の複数の画像をレジスタする方法。
【請求項２０】
前記形状保存特性は、基礎になる真の画像値を最小限に変更しながら雑音を除去する、請求項９から１９のいずれか１項に記載の複数の画像をレジスタする方法。
【請求項２１】
前記関数は、微分器をさらに含むことができる、請求項９から２０のいずれか１項に記載の複数の画像をレジスタする方法。
【請求項２２】
前記周波数領域内で、前記微分器の周波数応答は、高周波数の雑音を減衰させ、前記真の画像値を保存する、請求項９から２１のいずれか１項に記載の複数の画像をレジスタする方法。
【請求項２３】
前記微分器は、Ｓａｖｉｔｚｋｙ−Ｇｏｌａｙ微分器をさらに含むことができる、請求項９から２２のいずれか１項に記載の複数の画像をレジスタする方法。
【請求項２４】
整数シフトは、２つの画像を１／２画素以内にレジスタするのに使用される、請求項９から２３のいずれか１項に記載の複数の画像をレジスタする方法。
【請求項２５】
前記関数の前記周波数領域表現は、事前に計算される、請求項９から２４のいずれか１項に記載の複数の画像をレジスタする方法。
【請求項２６】
複数の画像内の前記整数部分と前記副画素部分との両方でのシフトの前記計算は、ハードウェアで実施される、請求項９から２５のいずれか１項に記載の複数の画像をレジスタする方法。
【請求項２７】
前記複数の画像の間で、前記方法は、
ａ）特徴抽出特性を含む前記関数を適用するステップと、
ｂ）各画像内の複数の点で画像特性を入手するステップと
をさらに含む、請求項９から２６のいずれか１項に記載の複数の画像をレジスタする方法。
【請求項２８】
前記特徴抽出特性は、雑音除去特性および／または平滑化特性をさらに含むことができる、請求項９から２７のいずれか１項に記載の複数の画像をレジスタする方法。
【請求項２９】
前記特徴抽出特性関数は、ハードウェアで実施される、請求項９から２８のいずれか１項に記載の複数の画像をレジスタする方法。
【請求項３０】
ａ）複数の画像を受け取るための画像入力手段と、
ｂ）レジスタされた画像を表示しまたは記憶するための出力手段と、
ｃ）前記複数の画像をレジスタするためのプロセッサと
を含み、前記プロセッサは、先行する請求項１〜２９のうちの任意の１つまたは複数に記載の前記方法を適用する
画像レジストレーション装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、画像レジストレーション方法に関し、具体的には、任意の次元の画像または信号の間の並進シフトを識別する画像レジストレーション方法に関する。
【背景技術】
【０００２】
複数の方法が、２次元（２Ｄ）画像および３次元（３Ｄ）画像をレジスタするために提案されてきた。画像レジストレーション技法は、様々な領域で幅広く使用されている。副画素画像レジストレーションおよび副画素変形測定は、動き推定および追跡、画像アライメント、医療画像処理、超解像再構成、画像モザイク化、衛星画像分析、および変化検出など、画像の正確なレジストレーションまたは小さい変形の正確な測定が要求される応用例で特に関心を持たれている。副画素画像レジストレーション技法は、産業応用例および医療応用例で表面検査、変形測定、および歪み測定にも使用される。
【０００３】
一般に、画像レジストレーション技法は、シーンまたは物体の複数の画像が異なる時におよび／または異なる位置から獲得される時に必要になる。課題は、画像が同一の座標系内にあるように画像を変換することに対する。これは、異なる画像からのデータの間の比較またはそのデータの統合を可能にする。画像の間の差は、並進、回転、スケーリング、または他のより複雑な変換を含む変換とすることができる。
【０００４】
２つの主要な手法が、副画素画像レジストレーションに関して使用可能である。第１の手法では、整数シフトおよび副画素／サブボクセル・シフトが、２つの別々のステップで見つけられる（すなわち、２ステップ法）。第２の主要な手法は、レジストレーションをコスト関数として扱い、副画素シフトは、ニュートン−ラフソン（ＮＲ）法および勾配ベースの（ＧＢ）方法などの方法に基づく反復誤差最小化プロセスを使用して見つけられる。画像を副画素精度でレジスタする方法の多くが、２ステップ法である。整数シフトを見つける最も一般的な方法は、２つの画像の相互相関（ＣＣ）または正規化相互相関（ＮＣＣ）内のピークを見つけることである。副画素／サブボクセル・シフトは、空間領域または周波数領域でのアップサンプリング、関数へのあてはめ、位相ベースの方法の使用、ならびに反復法および形状関数の使用など、様々な方法によって見つけられる。
【０００５】
デジタル画像相関（ＤＩＣ）は、画像内の副画素変形を見つけるための、最もよく知られた２ステップ法である。ＤＩＣは、異方性弾性膜の２Ｄ表面変形測定、人間の腱組織内のひずみ場の発見、回転ブレードの３Ｄ運動の測定、および金属シート溶接中の３Ｄ変形定量化を含む多数の応用例を有する。ＤＩＣ技法は、２Ｄおよび３Ｄの表面変形と体積変形とを測定することができる。以下のセクションでは、これらの技法を説明する。
【０００６】
２Ｄおよび３ＤのＤＩＣ
現在の２Ｄ／３Ｄ  ＤＩＣ技法は、複数の問題を有する。たとえば、整数シフトを見つけるという問題を検討する場合に、相互相関（ＣＣ）は十分に正確ではなく、正規化相互相関（ＮＣＣ）は計算集中型である。位相相関（ＰＣ）は、フーリエ・シフト・プロパティおよび正規化されたクロスパワー・スペクトルに基づいて整数シフトを見つける別の方法である。ＰＣは、雑音および強度の均一な変化に対してより頑健であり、ＣＣよりよく機能することが示されたが、十分な精度または計算効率を有しない。勾配相関（ＧＣ）法および正規化ＧＣ（ＮＧＣ）法は、ＣＣ関数の強度値の代わりに、画素ごとの中心差分近似に基づいて、複素値（実数値および虚数値）を使用する。
【０００７】
２つの画像の間の副画素シフトは、複数の方法によって見出され得る。たとえば、空間領域での補間が、時々使用されるが、これは、通常は計算集中型である。この問題に対処するために、周波数領域での補間（ＦＦＴアップサンプリング）が提案された。しかし、この方法の精度は、アップサンプリング比ならびに補間方法によって制限され、高いサンプリング比に関して低速である。曲線あてはめは、副画素シフトを見つけるためにＣＣ関数またはＧＣ関数のピーク付近で関数を数値的にあてはめるのに使用され得る。関数の例は、ガウス関数、二次関数、および修正メキシカン・ハット・ウェーブレット（ｍｏｄｉｆｉｅｄ  Ｍｅｘｉｃａｎ  ｈａｔ  ｗａｖｅｌｅｔ）関数を含む。しかし、これらの方法の精度は、ＣＣ関数形状またはＧＣ関数形状に依存する。さらなる欠点は、これらの方法が、通常、最良のあてはめを見つけるための、時間のかかる反復最適化プロセスを含むことである。
【０００８】
別の技法は、フーリエ周波数シフト理論および２つの画像の強度値の位相差の傾きに基づく、位相ベースの手法である。この方法では、位相アンラッピングが、１画素より大きいシフトについて必要である。以前の技法は、まず、整数シフトを判定し、その後、位相差データに対して２Ｄ回帰を使用して、副画素シフトを見つけた。この方法は、エイリアシングに敏感であり、画像が、第１ステップで１画素の半分以内にレジスタされることを必要とする。ＤＩＣおよび実験力学のさらなる技法は、副画素変位を測定するためのＮＲおおびＧＢなどの反復法を使用する。これらは、高い計算コストを有し、したがって低速である。これらは、強度値の差に基づくので、これらの方法は、強度変化に敏感であり、小さいシフトのみに適するものになる。これは、多くの応用例にとって重大な制限である。さらに、反復ＮＲ法の各サブセット変形は、補間を伴い、これが、追加の系統誤差を導入する。
【０００９】
体積変形
デジタル体積相関（ＤＶＣ）（体積ＤＩＣとも称する）は、体積における３Ｄサブボクセル変形（または変位）を測定するためのＤＩＣの拡張である。ＤＶＣは、計算モデルに情報を与えるために、脳などの人間の軟組織の変位を定量化するのに使用されてきた。さらに、ＤＶＣは、材料の構造的挙動および機械的挙動を査定する、幅広く使用されている技法の本質的部分であり、この技法では、デバイスが変形前および変形中の画像を記録する間に、試験材料が、印加された力の下で変形される。Ｘ線コンピュータ断層撮影（ＣＴ）は、体積変形を取り込むのに使用される最も一般的な画像診断法である。たとえば、ＤＶＣ技法が、フォーム［１３］、合成物、骨、骨格インプラント、ｐｏｌｙｍｅｒ  ｂｏｎｄｅｄ  ｓｕｇａｒ、およびパンくずの機械的挙動を査定するためにＸ線ＣＴまたはマイクロＣＴに適用された。ＤＶＣ技法は、Ｘ線ＣＴ画像および放射光ラミノグラフィ画像での亀裂の開始および伝搬の調査ならびにＸ線ＣＴ画像での熱特性の測定にも使用された。しかし、低品質３Ｄイメージング技術でのＤＶＣの使用は、現在のＤＶＣ技術の限界ゆえに制限される。たとえば、ＤＶＣは、光干渉断層撮影（ＯＣＴ）画像から弾性スティフネスを測定するのに、２０１３年に初めて使用された。
【００１０】
ＳｕｔｔｏｎおよびＨｉｌｄは、最近の論文で、正確なグレイ・レベル補間の実行、適切な形状関数の選択、および短時間での膨大な量のデータの処理を含む現在のＤＶＣ技法の課題の一部を要約した。さらに、ＤＶＣアルゴリズムは、パラメータのよい初期推定値を必要とする。これらの課題は、ＤＶＣに使用される既存技法固有の限界から生じる。既存ＤＶＣ技法のほとんどは、整数シフトを見つけるのに体積の３Ｄ−ＣＣを、サブボクセル・シフトを見つけるのに反復非線形最適化を使用する。ＤＶＣアルゴリズムの計算コストは、Ｇ画素（ボクセル）の格子点間隔および画像にわたる（体積の内部の）Ｓ画素２（ボクセル３）のサブセット・サイズに関して２Ｄ−ＤＩＣアルゴリズムよりＧ×Ｓ倍高く、したがって低速である。たとえば、以前に、８プロセッサを有する並列コンピューティング・アーキテクチャは、（３９×３９×３９）点からなる格子内の、サイズ４１ボクセル×４１ボクセル×４１ボクセルのＤＶＣを計算する計算時間を、４５．７時間から５．７時間に短縮することしかできなかった。この限界は、高解像度画像および大量のデータに関してより問題になる。さらに、ＤＶＣアルゴリズムで使用されるＣＣは、雑音および画像照明の変化に敏感であり、貧弱なテクスチャを有するか大きい変形を受けた画像に関して失敗する。したがって、ＤＶＣの使用は、豊かなテクスチャを有する体積内の小さい変形の測定に制限された。レジストレーションを改善するために、ＣＣは、いくつかの方法でＮＣＣに置換された。しかし、ＮＣＣは、画像照明の変化を扱う際にＣＣを越えるいくつかの利点を有するが、ＣＣの複数の他の限界には対処しない。たとえば、ＮＣＣは、実質的にＣＣより計算的に要求が厳しく、大きい変形に関して貧弱に動作する。Ｐａｎ他（Ｐａｎ，  Ｂ．，  Ｙｕ，  Ｌ．，  Ｗｕ，  Ｄ．  Ｈｉｇｈ−ａｃｃｕｒａｃｙ  ２Ｄ  ｄｉｇｉｔａｌ  ｉｍａｇｅ  ｃｏｒｒｅｌａｔｉｏｎ  ｍｅａｓｕｒｅｍｅｎｔｓ  ｕｓｉｎｇ  ｌｏｗ−ｃｏｓｔ  ｉｍａｇｉｎｇ  ｌｅｎｓｅｓ：  ｉｍｐｌｅｍｅｎｔａｔｉｏｎ  ｏｆ  ａ  ｇｅｎｅｒａｌｉｚｅｄ  ｃｏｍｐｅｎｓａｔｉｏｎ  ｍｅｔｈｏｄ．  Ｍｅａｓｕｒｅｍｅｎｔ  Ｓｃｉｅｎｃｅ  ａｎｄ  Ｔｅｃｈｎｏｌｏｇｙ，  （２０１４年），  ２５：２）は、大きい変形を測定するために基準画像を更新する増分ＤＩＣ法を提案することによって、後者の限界に対処したが、彼らの方法は、計算時間をかなり増加させた。
【００１１】
ＣＣおよびＮＣＣの制限は、勾配相関（ＧＣ）および位相相関（ＰＣ）など、いくつかの２Ｄ法によって対処された。ＧＣは、２つの座標方向の強度値の中心差分を組み合わせて、一方の実数部分画像
【数１】
2019517073.tif 000003
を乗算することと、これを他方の実数部分画像に加算することとによって単一の複素画像を形成する。この手法は、２つの実数値の情報を単一の複素数値として符号化することを可能にする。ＰＣは、２つの画像の強度の正規化されたクロスパワー・スペクトルを使用して、その相対シフトを見つける。ＰＣおよびＧＣは、画像の強度値を直接には使用しないので、両方とも、貧弱なテクスチャを有する画像内でシフトを見つけることにおいて、ＣＣより頑健である。ＧＣは、後に２Ｄ副画素レジストレーション・アルゴリズム内でも使用された。
【００１２】
いくつかの手法が、試験体積が大きい変形を有する場合にＤＶＣ応用例でのＣＣおよびＮＣＣの短所に対処するために提案された。ストレッチ相関（ｓｔｒｅｔｃｈ−ｃｏｒｒｅｌａｔｉｏｎ）法が、大きい変形を測定する際のＤＶＣの制限に対処するために過去に提案された。ストレッチ相関は、体積の高速フーリエ変換（ＦＦＴ）を使用して、フーリエ領域で実施された。ストレッチ相関は、小さい回転およびせん断を仮定して、変形勾配テンソルを直交回転および対称右ストレッチ・テンソルに分解することによって、部分画像のストレッチを極座標系で考慮に入れる。しかし、ストレッチ相関法は、部分画像が変位されまたはシフトされる時ではなく、体積がストレッチされる場合の大きい変形でのレジストレーション性能を改善することだけができる。最近、Ｂａｒ−Ｋｏｃｈｂａ他（Ｂａｒ−Ｋｏｃｈｂａ，  Ｅ．，  Ｔｏｙｊａｎｏｖａ，  Ｊ．，  Ａｎｄｒｅｗｓ，  Ｅ．  Ａ  Ｆａｓｔ  Ｉｔｅｒａｔｉｖｅ  Ｄｉｇｉｔａｌ  Ｖｏｌｕｍｅ  Ｃｏｒｒｅｌａｔｉｏｎ  Ａｌｇｏｒｉｔｈｍ  ｆｏｒ  Ｌａｒｇｅ  Ｄｅｆｏｒｍａｔｉｏｎｓ．  Ｅｘｐｅｒｉｍｅｎｔａｌ  Ｍｅｃｈａｎｉｃｓ，  ２０１５年，  ５５：  ２６１）は、大きい変形を測定するための反復ＤＶＣ手法およびＣＣ係数の重み付け関数を提案した。Ｂａｒ−Ｋｏｃｈｂａ他の方法は、２Ｄ応用例に関してＰａｎ他によって提案されたものに似た手法を使用し、この手法は、大きい部分画像から開始し、変形を測定し、２つの体積を歪曲させ、体積の間の誤差を測定し、より小さい部分画像を用いて別の反復を継続すべきか、またはプロセスを停止すべきかのいずれであるのかを誤差しきい値に基づいて判断するというものであった。Ｂａｒ−Ｋｏｃｈｂａ他の方法での重み付け関数の使用の目的は、高周波数に重みを付けることによって変位場の分解能を高めることであった。Ｂａｒ−Ｋｏｃｈｂａ他の方法の重み付け関数は、粒子画像速度測定に関して提案されたＮｏｇｕｅｉｒａ他の方法（Ｎｏｇｕｅｉｒａ，  Ｊ．，  Ｌｅｃｕｏｎａ，  Ａ．，  Ｒｕｉｚ−Ｒｉｖａｓ，  Ｕ．，  Ｒｏｄｒｉｇｕｅｚ，  Ａ．  Ａｎａｌｙｓｉｓ  ａｎｄ  ａｌｔｅｒｎａｔｉｖｅｓ  ｉｎ  ｔｗｏ−ｄｉｍｅｎｓｉｏｎａｌ  ｍｕｌｔｉｇｒｉｄ  ｐａｒｔｉｃｌｅ  ｉｍａｇｅ  ｖｅｌｏｃｉｍｅｔｒｙ  ｍｅｔｈｏｄｓ：  ａｐｐｌｉｃａｔｉｏｎ  ｏｆ  ａ  ｄｅｄｉｃａｔｅｄ  ｗｅｉｇｈｔｉｎｇ  ｆｕｎｃｔｉｏｎ  ａｎｄ  ｓｙｍｍｅｔｒｉｃ  ｄｉｒｅｃｔ  ｃｏｒｒｅｌａｔｉｏｎ．  Ｍｅａｓｕｒｅｍｅｎｔ  Ｓｃｉｅｎｃｅ  ａｎｄ  Ｔｅｃｈｎｏｌｏｇｙ，  ２００２年，  １３  ９６３−７４）から適合された。各反復でＣＣ出力から異常値を識別し、除去し、最終的な相互相関出力のピークに二変数ガウス関数をあてはめることによって彼らの方法でサブボクセル・シフトを見つけた。Ｂａｒ−Ｋｏｃｈｂａ他［３４］の反復法での重み付け関数の使用は、大きいシフトでのＣＣの性能を改善するが、ＣＣの低域通過挙動を除去することはできない。さらに、Ｂａｒ−Ｋｏｃｈｂａ他［３４］の手法は、計算集中型の反復プロセスである。
【００１３】
体積勾配相関を使用する別の方法が、低コントラスト画像と２Ｄ投影された体積との間で２Ｄ−３Ｄレジストレーションを実行する時のＣＣの短所を克服するために提案された。勾配相関は、２Ｄ投影された体積の勾配画像の２つの座標のＮＣＣ値の平均値を使用する。体積勾配相関は、医療画像での２Ｄ−３Ｄレジストレーションに関してＮＣＣよりよく機能することが示され、ＣＣが失敗する、臨床３Ｄ  ＣＴデータを２Ｄ  Ｘ線画像にレジスタすることができる。体積勾配相関は、低コントラスト応用例ではＣＣおよびＮＣＣよりよく機能するが、２Ｄの事例にのみ適用可能であり、ＤＶＣには使用され得ない。しかし、低コントラスト体積の処理は、ＤＶＣ応用例での大きな課題である。というのは、２Ｄの事例で実行されるように、ランダム・スペックル・パターンを追加することによってコントラストを高めることが、難しいからである。
【００１４】
したがって、現在の技術は、低テクスチャ画像に含まれる、多数の２Ｄ／３Ｄ応用例に関する適切なレベルの精度、速度、雑音に対する頑健性、および／または計算効率を提供しない。特定の問題は、諸方法の精度が、多くの応用例で非常に正確な副画素レジストレーションを可能にするのに十分ではないことである。
【発明の概要】
【発明が解決しようとする課題】
【００１５】
本発明の目的は、従来のシステムの効率、および／または精度、および／または頑健性を向上させる画像レジストレーション方法を提供することである。この画像レジストレーション方法は、任意の次元の画像に適用可能とすることができる。
【００１６】
本発明の目的は、少なくとも既存システムの不利益をある程度克服でき、または少なくとも既存システムに対する有用な代替案を提供する、画像レジストレーション方法を提供することである。
【００１７】
本発明のさらなる目的は、少なくとも既存システムの不利益を程度克服できる、副画素精度を有する画像レジストレーション方法を提供することである。
【００１８】
本発明のさらなる目的は、以下の説明から明白になる。
【課題を解決するための手段】
【００１９】
したがって、第１の態様では、本発明は、任意の次元の複数の画像の間での画像レジストレーション方法に存するとおおまかに言うことができる。この方法は、既存方法と比較して、整数部分と副画素部分との両方でシフトを見つけることの精度、速度、および頑健性を改善する２ステップ方法である。整数部分では、フィルタが、画像情報を強調し、雑音の影響を低下させるために画像に適用される。副画素部分では、副画素シフトを見つける位相ベースの方法に含まれる雑音およびスプリアス周波数成分の影響を低下させるために、複数の方法が提案された。本方法は、まず、幅広い形で説明され、その後、特定の２Ｄ実施態様が説明される。
【００２０】
幅広い形の画像レジストレーション方法
本発明は、任意の次元の複数の画像の間での画像レジストレーション方法であって、
１つまたは複数のフィルタ数の周波数領域表現を入手することであって、各フィルタ関数は、少なくとも１つの画像特性を強調する、入手するステップ
を含む方法に存するとおおまかに言うことができる。特性は、低減された雑音または改善された信号対雑音比とすることができる。
【００２１】
一実施形態では、フィルタ関数は、カーネルとして空間領域で定義される。一実施形態では、本方法は、フィルタ関数のフーリエ変換をとるステップを含む。
【００２２】
一実施形態では、単一の関数を形成するために、フィルタ関数が組み合わされる。一実施形態では、組合せは、周波数領域での乗算のステップを含む。
【００２３】
一実施形態では、本方法は、フィルタ関数を二乗するステップを含む。
【００２４】
一実施形態では、本方法は、複数の画像のフーリエ領域相互相関を入手するステップを含む。
【００２５】
一実施形態では、本方法は、二乗されたフィルタ関数と複数の画像のフーリエ領域相互相関とを組み合わせることによって、フィルタリングされた相互相関（ＦＣＣ）を形成するステップを含む。
【００２６】
一実施形態では、画像は、１つ、２つ、３つ、またはそれより高い次元を有する。
【００２７】
一実施形態では、ＦＣＣで使用される１つまたは複数のフィルタ関数は、複数の次元を有する。好ましくは、フィルタ関数は、画像と同一個数の次元を有する。
【００２８】
一実施形態では、本方法は、ＦＦＣ値に重みを付けるステップを含む。
【００２９】
一実施形態では、本方法は、フーリエ領域のフィルタリングされた相互相関の大きさを入手するステップを含む。一実施形態では、大きさは、最大ピークを見つけるのに使用される。
【００３０】
一実施形態では、整数シフトは、２つの画像を１画素未満以内でレジスタするのに使用される
【００３１】
一実施形態では、副画素シフトは、以前のステップで見つけられた整数シフト値を使用してレジスタされた画像の位相差データを使用して計算される。
【００３２】
一実施形態では、ＦＣＣで使用されるフィルタ関数は、強度勾配を抽出する微分カーネルと、雑音の影響を低減する平滑化関数とを含む。
【００３３】
一実施形態では、フィルタ関数は、フーリエ領域で適用される。一実施形態では、本システムは、複数の異なるカーネル／関数を含む。
【００３４】
好ましくは、ＦＣＣで使用されるフィルタ関数の周波数領域表現は、事前に計算される。
【００３５】
したがって、さらなる態様によれば、本発明は、複数の画像の間での画像レジストレーション方法であって、
関数（フィルタ）の周波数領域表現を入手するステップと、
関数（フィルタ）の周波数領域表現を相互相関の周波数領域表現に適用するステップと
を含む方法に存するとおおまかに言うことができる。
【００３６】
一実施形態では、周波数領域表現を入手するステップは、時間領域関数（フィルタ）または空間領域関数（フィルタ）からフーリエ変換をとることを含む。
【００３７】
好ましくは、関数は、複数の可能なフィルタ関数のうちの１つである。
【００３８】
好ましくは、本方法は、複数の事前に計算された関数から関数を選択するステップを含む。
【００３９】
好ましくは、本方法は、第２のまたは複数の異なる関数を用いて本方法を繰り返すステップを含む。
【００４０】
好ましくは、関数は、平滑化関数および／または微分関数をさらに含むことができる。
【００４１】
好ましくは、平滑化関数は、平滑化カーネルを含む。好ましくは、平滑化関数は、雑音の影響を低減させる。好ましくは、平滑化関数は、少なくとも１つの形状保存特性を含む。
【００４２】
好ましくは、本方法は、複数の画像に適用され、事前に計算された平滑化関数が、再利用される。
【００４３】
好ましくは、周波数領域表現は、微分演算をさらに含む。好ましくは、微分演算または微分演算子は、画像特徴のうちの少なくとも１つを強調し、かつ／または抽出する。
【００４４】
好ましくは、平滑化フィルタの周波数領域表現は、相互相関に適用される前に変更される。好ましくは、変更は、画像内容に依存する。
【００４５】
好ましくは、平滑化フィルタの周波数領域表現は、重み付けプロファイルを含む。一実施形態では、本方法は、周波数領域表現に重みを付けるステップを含む。
【００４６】
好ましくは、重み付けプロファイルは、二乗された値である。一実施形態では、フィルタリングされた相互相関の出力は、総フィルタリングされた相互相関（ＦＣＣ）を計算する前に重みを付けられる。
【００４７】
好ましくは、本方法は、より滑らかな関数の周波数領域表現を相互相関の表現に適用することによって、平滑化された相互相関の周波数領域表現を入手するステップを含む。
【００４８】
好ましくは、適用は、表現の乗算を含む。
【００４９】
好ましくは、相互相関は、少なくとも部分的に、第１の画像および第２の画像の複素共役の乗算によって計算される複数の画像の相互相関を含む。
【００５０】
好ましくは、本方法は、複数の画像のうちの少なくとも１つに窓関数を適用するステップを含む。
【００５１】
好ましくは、周波数領域表現のうちの少なくとも１つは、フーリエ表現である。
【００５２】
好ましくは、複数の画像は、２つ以上の部分画像を含む。
【００５３】
好ましくは、本方法は、複数の画像の表現に周波数領域変換を適用するステップを含む。
【００５４】
好ましくは、本方法は、平滑化された相互相関の周波数領域表現の空間領域表現を計算するステップを含む。
【００５５】
好ましくは、本方法は、ベクトル値を有する画像として複数の空間領域表現を組み合わせるステップを含む。
【００５６】
好ましくは、本方法は、２つの画像の間の相対変位を推定するステップを含む。
【００５７】
したがって、さらなる態様では、本発明は、複数の画像の間での画像レジストレーション方法であって、整数シフトを見つけるために周波数領域の複数の画像の間のフィルタリングされた相互相関（ＦＣＣ）を入手するステップを含む方法に存するとおおまかに言うことができる。
【００５８】
一実施形態では、画像レジストレーション方法は、平滑化微分器カーネルなど、画像特徴を抽出しもしくは強調し、かつ／または雑音の影響を低減するフィルタを適用するステップを含む。
【００５９】
好ましくは、フィルタリングされた相互相関（ＦＣＣ）を入手するステップは、複数の画像に対する相互相関によって１つまたは複数のフィルタ関数を乗算することを含む。
【００６０】
好ましくは、１つまたは複数のフィルタ関数は、事前に計算される。  
【００６１】
好ましくは、フィルタ関数は、平滑化微分器関数である。
【００６２】
好ましくは、本方法は、時間領域または空間領域において、フィルタリングされた相互相関表現を入手するステップを含む。
【００６３】
好ましくは、本方法は、周波数領域のフィルタリングされた相互相関（ＦＣＣ）の出力から画像の間の相対整数変位の推定値を入手するステップを含む。
【００６４】
好ましくは、雑音およびエイリアシングの影響は、副画素部分内の位相データから除去された。
【００６５】
さらなる態様によれば、本発明は、複数の画像の間での画像レジストレーション方法であって、
各画像内の複数の点での画像特性を入手するステップと、
画像にフィルタ／カーネル（平滑化勾配フィルタを含むことができる）を適用するステップと
を含む方法に存するとおおまかに言うことができる。
【００６６】
したがって、さらなる態様では、本発明は、
複数の画像を入手するイメージャと、
レジスタされた画像を表示しまたは記憶する出力と、
複数の画像をレジスタするプロセッサと
を含み、プロセッサは、上記の態様または実施形態のうちの任意の１つまたは複数に記載の方法を適用する、画像レジストレーション装置に存するとおおまかに言うことができる。
【００６７】
本発明は、添付の特許請求の範囲のいずれか一項に記載の画像レジストレーション方法または画像レジストレーション装置にも存することができる。
【００６８】
２Ｄ画像レジストレーションのための画像レジストレーション方法の実施形態
本発明は、より正確な２Ｄ画素並進シフトを提供する。これは、勾配測定の複数の隣接する点を組み合わせる、勾配の計算によって達成される。このプロセスは、平滑化関数などの特徴抽出関数の適用をも含むことができる。これは、微分器カーネル（勾配）など、画像特性の一部を抽出しまたは強調するさらなる演算子と組み合わされ（またはこれを含み）得る。勾配は、計算における最小限の要因と思われ、より複雑なまたはより高次の関数は、計算負荷を増加させるが、平滑な勾配または微分器カーネルと組み合わされた平滑化フィルタの追加は、画像レジストレーションに対する大きい有益な影響を有する。
【００６９】
一実施形態では、勾配推定は、演算子をさらに含む。一実施形態では、演算は、画像特性のうちの少なくとも１つを強調する。
【００７０】
一実施形態では、特徴抽出関数は、平滑化関数を含む。
【００７１】
一実施形態では、画像特性は、強度を含む。
【００７２】
一実施形態では、画像特性を入手するステップは、抽出を介するものとすることができる。
【００７３】
一実施形態では、画像特性のセットが、入手され得る。
【００７４】
一実施形態では、特性は、画像の複数の次元またはすべての次元で抽出され得る。
【００７５】
一実施形態では、特徴抽出関数は、少なくとも１つの形状保存特性を含む。
【００７６】
一実施形態では、形状保存特性は、基礎になる真の値を最小限に変更しながら雑音を除去する。
【００７７】
一実施形態では、特徴抽出関数は、少なくとも１つの雑音除去特性を含む。
【００７８】
一実施形態では、勾配推定は、雑音低減プロパティを有する。
【００７９】
一実施形態では、平滑化関数は、多項式への勾配のあてはめに基づく。
【００８０】
一実施形態では、特徴抽出関数は、周波数領域で画像に適用される。
【００８１】
一実施形態では、勾配は、複数の勾配値を含む。
【００８２】
一実施形態では、あてはめる手段は、移動最小二乗（ｒｕｎｎｉｎｇ  ｌｅａｓｔ−ｓｑｕａｒｅ）である。
【００８３】
一実施形態では、平滑化関数の周波数応答は、要求される形状を有する。
【００８４】
一実施形態では、特徴抽出器関数は、要求される形状を有する。
【００８５】
一実施形態では、平滑化関数および／または特徴抽出器関数は、畳み込みカーネルを使用して実施される。一実施形態では、特徴抽出器関数は、フィルタである。
【００８６】
一実施形態では、平滑化関数および／または特徴抽出器関数は、ハードウェアで実施される。
【００８７】
一実施形態では、勾配推定は、Ｓａｖｉｔｚｋｙ−Ｇｏｌａｙ微分器を使用する。
【００８８】
一実施形態では、複数の点が、画像の画素を表す。
【００８９】
一実施形態では、平滑化関数は、画像または画像表現の行および列に適用される。一実施形態では、平滑化関数は、画像または画像表現の複数の次元またはすべての次元に適用される。
【００９０】
一実施形態では、勾配は、画像の複数の次元またはすべての次元で推定される。一実施形態では、次元は、直交方向である。
【００９１】
一実施形態では、特徴抽出器関数は、画像の複数の次元またはすべての次元に適用される。一実施形態では、これは、周波数領域にある。
【００９２】
一実施形態では、本方法は、複数の画像を入手するステップを含む。
【００９３】
一実施形態では、画像は、カメラまたはビデオ・カメラによって入手される。一実施形態では、画像は、イメージング・システム（たとえば、ＭＲＩ、ＣＴスキャン、衛星、カメラまたはビデオ・カメラ）から入手される。
【００９４】
一実施形態では、本方法は、画像を周波数領域に変換するステップを含む。好ましくは、変換は、空間領域からの変換である。好ましくは、フーリエ変換が使用される。一実施形態では、画像は、それらの間のシフトを見つけるために、周波数領域に変換され、操作され、周波数領域から逆変換される。
【００９５】
一実施形態では、本方法は、相互相関を計算するステップを含む。一実施形態では、相互相関は、正規化される。
【００９６】
一実施形態では、本方法は、相互相関を計算した。
【００９７】
一実施形態では、本方法は、窓関数を適用するステップを含む。好ましくは、窓関数は、画像の高周波数情報を保存する。
【００９８】
一実施形態では、本方法は、整数画素シフトを計算する。
【００９９】
さらなる態様では、本発明は、複数の画像の間の画像レジストレーション方法であって、
各画像内の複数の点での画像特性を入手するステップと、
画像特性の勾配を推定するステップであって、勾配推定は、形状保存関数を含む、推定するステップと
を含む方法に存するとおおまかに言うことができる。
【０１００】
本発明は、より正確な画素並進シフトを提供する。これは、勾配測定の複数の隣接する点を組み合わせる、勾配の計算によって達成される。勾配は、計算における最小限の要因と思われ、より高次の関数は、計算負荷を増加させるが、形状保存勾配推定の追加は、正確な勾配推定値と雑音の除去とを作成する。好ましくは、形状保存関数は、雑音を除去し、基礎になる真の値を最小限に変更する。
【０１０１】
さらなる態様では、本発明は、複数の画像の間の画像レジストレーション方法であって、
各画像内の複数の点での画像特性を入手するステップと、
画像特性の雑音特性を入手するステップと、
雑音特性に基づいてサンプル・サイズを選択するステップと、
サンプル・サイズ内の画像特性データに基づいて画像特性の勾配を計算するステップと
を含む方法に存するとおおまかに言うことができる。
【０１０２】
一実施形態では、副画素シフト推定の画像特性は、位相差である。
【０１０３】
一実施形態では、位相差は、２つの画像の間の位相差である。
【０１０４】
一実施形態では、雑音特性は、標準偏差である。
【０１０５】
一実施形態では、雑音消去は、位相差データの標準偏差である。
【０１０６】
一実施形態では、本方法は、サンプル・サイズ以内の内の画像特性をフィルタリングするステップを含む。
【０１０７】
一実施形態では、フィルタリングするステップは、高周波数雑音を除去する。
【０１０８】
一実施形態では、方法は、特徴サイズ以内の内の画像特性をフィルタリングするステップを含む。
【０１０９】
一実施形態では、サンプル・サイズ選択は、雑音特性を含む式に基づいて適合される。
【０１１０】
一実施形態では、本方法は、位相差データを平滑化するために位相差をフィルタリングするステップを含む。
【０１１１】
一実施形態では、本方法は、スプリアス周波数を除去するために位相差をフィルタリングするステップを含む。
【０１１２】
一実施形態では、フィルタは、メディアン・フィルタである。一実施形態では、フィルタは、２Ｄメディアン・フィルタである。
【０１１３】
一実施形態では、本方法は、副画素シフトを判定する。
【０１１４】
一実施形態では、第２の態様は、第１の態様と組み合わされる。
【０１１５】
一実施形態では、サンプル・サイズは、位相データ・サイズである。
【０１１６】
さらなる態様では、本発明は、複数の画像の間の画像レジストレーション方法であって、
画像の間の整数画素シフトの推定値を入手するステップと、
画像の間の副画素シフトの推定値を入手するステップと
を含み、副画素シフトは、１画素未満だけ離れた画像を使用して計算される、方法に存するとおおまかに言うことができる。
【０１１７】
１画素未満だけ離れた画像の仮定または使用は、副画素計算の複雑さを低減させる。というのは、たとえば、画像位相が、処理の前にアンラッピングされる必要がないからである。これは、副画素シフトの計算に関するより清浄でより明瞭なデータをもたらす。
【０１１８】
一実施形態では、副画素シフトは、まず位相データをアンラッピングすることなく計算される。
【０１１９】
一実施形態では、画像は、１画素未満だけ離れていると仮定される。
【０１２０】
一実施形態では、画像の間の副画素シフトの推定値を入手するステップは、入手された整数画素シフトによって画像のうちの１つをシフトするステップを含む。
【０１２１】
一実施形態では、整数画素シフトおよび副画素シフトは、本明細書で説明される実施形態のうちの任意の１つである。
【０１２２】
さらなる態様では、本発明は、複数の画像の間の画像レジストレーション方法であって、
オフセット値を減算し、最大値によって除算することによって、複数の画像特性を正規化するステップ
を含む方法に存するとおおまかに言うことができる。
【０１２３】
好ましくは、オフセット値は、画像特性の平均値である。
【０１２４】
好ましくは、平均値は、平均ＤＣ値である。
【０１２５】
好ましくは、最大値は、画像特性の最大値である。
【０１２６】
好ましくは、上で説明された特徴は、可能な時に、説明された諸態様のいずれにも適用されると考えられなければならない。
【０１２７】
開示される主題は、部分、要素、または特徴のうちの２つ以上の任意のまたはすべての組合せで、個別にまたは集合的に本明細書で参照されまたは示された部分、要素、または特徴に存するとおおまかに言われ得る任意の次元の画像レジストレーション方法をも提供する。本発明が関係する技術分野で既知の同等物を有する特定の整数が、本明細書で述べられる場合に、そのような既知の同等物は、本明細書に組み込まれていると考えられる。
【０１２８】
そのすべての新規の態様において考慮されなければならない本発明のさらなる態様は、以下の説明から明白になる。
【０１２９】
本発明の複数の実施形態が、これから図面を参照して例によって説明される。
【図面の簡単な説明】
【０１３０】
【図１】本発明の２Ｄ形態での２ステージ方法を使用して画素シフトを識別するのに相互相関が使用される実施形態を示す流れ図である。
【図２】本発明の２Ｄ形態で整数画素シフトを識別するのにＳＧＤが使用される実施形態を示す流れ図である。
【図３】本発明の２Ｄ形態で副画素シフトを識別するのに使用される位相差詳細の２Ｄ回帰詳細が場合の実施形態を示す流れ図である。
【図４】本発明の２Ｄ形態での中心差分および７点三乗一次導関数Ｓａｖｉｔｚｋｙ−Ｇｏｌａｙ補間の周波数応答特性を示すプロットである。
【図５】高周波数成分の微細画像詳細に対する雑音の分布を示す、（ａ）元の画像に関する、および（ｂ）少量の白色ガウス雑音を加えた、サンプル画像（図９ａ）のパワー・スペクトルを示す２Ｄヒート・マップである。
【図６】Ｈａｍｍｉｎｇ、Ｈａｎｎ、Ｂｌａｃｋｍａｎ、およびＴｕｋｅｙ（ａ＝０．２５）の窓関数の周波数応答を示す図である。
【図７】本発明の２Ｄ形態での正規化手順を示す図である。
【図８】画像レジストレーションに関する標準ランドサット画像を示す図である。
【図９】（ａ）パリのランサット画像上の部分画像の中心点。（ｂ）サンプル部分画像（１２８画素×１２８画素）を示す図である。
【図１０】（ａ）我々の方法および既存の方法に関するパリのランドサット画像（１２８画素×１２８画素）の４００個の部分画像内のシフトを推定する際の平均レジストレーション誤差（画素）の実施形態と、（ｂ）我々の方法の整数部分の誤差メトリックとして０．８５より大きい正規化されたＳＧＧＣ値（９）を有する点の平均個数の実施形態と、（ｃ）一実施形態の誤差メトリックとして位相差データに対する平面あてはめ誤差の実施形態とを示すプロットである。
【図１１】（ａ）パリのランドサット画像と（ｂ）ガウス雑音が追加されたパリのランドサット画像とを示す画像である。
【図１２】（ｂ）同一サイズの画像の１°回転の理想的回転パターンのヒート・マップを示す図である。回転パターンは、我々の方法（ｃ）、既存の方法によって見つけられ、（ａ）ミシシッピのランドサット画像の中心付近の回転に関して本発明の２Ｄ形態を使用して判定された変位ベクトル。
【図１３】本発明の２Ｄ形態の実施形態を示す流れ図である。
【図１４】本発明の２Ｄ形態で整数シフトを見つける、本発明の実施形態を示す流れ図である。
【図１５】整数シフトを見つける、本発明の広義の形態の実施形態を示す流れ図である。
【図１６】整数シフトを見つける、本発明の広義の形態の第２の実施形態を示す流れ図である。
【図１７】整数シフトを見つける、本発明の広義の形態の第３の実施形態を示す流れ図である。
【図１８】副画素シフトを見つける、本発明の広義の形態の実施形態を示す流れ図である。
【図１９】副画素シフトを見つける、本発明の広義の形態の第２の実施形態を示す流れ図である。
【図２０】副画素シフトを見つける、本発明の広義の形態の第３の実施形態を示す流れ図である。
【図２１】本方法の実施形態の効率を示すプロットである。
【図２２】画像レジストレーション方法を示す概略図である。
【発明を実施するための形態】
【０１３１】
この説明全体を通じて、同様の符号が、異なる実施形態内の同様の特徴を参照するのに使用される。
【０１３２】
例の実施形態では、画像レジストレーションは、医療分野で使用され得る。たとえば、本システムは、首の頸静脈の変動または膨張を評価するのに使用され得る。２つのカメラを使用する立体視システムが、接触なしで静脈を３Ｄで監視するのに使用され得る。立体視システムが説明されるが、さらなるカメラまたは画像ソースが使用され得る。提案される方法では、光画像が、患者の首に投影される。２つのカメラが使用されるので、３Ｄ画像が、これらの組合せを介して形成され得る。しかし、一方または両方のカメラから入手される画像が、ある時間期間にわたって画像の間のすべての差の比較のために正確に整列されることが重要である。たとえば、測定中にユーザの首の運動がある場合には、システムは、静脈運動が患者運動から区別され得るように、運動を推定できなければならない。
【０１３３】
より広範に、画像は、異なる時に撮影された同一タイプの画像内の差を探すために比較され得る。たとえば、造影剤（マーカを用いる血管造影の、トレーサを用いるＭＲＩに関する技法でなど）の使用の前後の位置のマッピングである。第２の例では、画像レジストレーションは、腫瘍成長または変性疾患など、構造変化をマッピングすることができる。
【０１３４】
さらなる例では、本発明の画像レジストレーションは、脳刺激を伴う機能的ＭＲＩまたはトレーサを用いるＰＥＴイメージングなど、機能の変化を検出するのに使用され得る。
【０１３５】
一般に、画像レジストレーションは、写真を撮影しまたは他の形で画像を入手しつつあり、患者（または他の物体）が定位置に適切に固定され得ないシステム内で使用され得る。
【０１３６】
さらなる例では、本発明の画像レジストレーションは、異なる画像ソースからの情報を関係付けるのに使用され得る。これは、画像を入手するのに使用される異なる装置を介するものとするか、同一のまたは異なる機械によって行われた異なる測定を関係付けることとすることができる。画像は、２次元、３次元、もしくは多次元で動作するイメージャまたは画像入力手段によって入手され得る。たとえば、イメージャ／画像入力手段は、カメラ、ＭＲＩ機械、またはｘ線機械とすることができる。
【０１３７】
図１５、図１６、および図１７は、方法のワークフローが、少なくとも、３（図１５）での平滑化関数１４２の別々の計算または（図１６）での任意の次元の任意の個数のフィルタおよび図１７での任意の次元の平滑化微分器カーネル）を可能にするように適合される、本発明の広義の形態の代替実施形態を示す。図１５の平滑化関数１４２、図１６のフィルタ（１７１）、および図１７の平滑化微分器カーネル（１７２）は、特定の画像もしくは部分画像の知識を獲得することによって選択され得、または、包括的な平滑化関数が選択され得る。
【０１３８】
周波数変換が、部分画像の次元と同等の周波数領域表現を作るために適用される（たとえば、ＤＦＴ）。その後、平滑化関数の周波数領域表現が、微分演算子ｉ（ｗｋ））１５３によって乗算されて、各方向ベクトル１４３内の平滑化関数の勾配に対する同等物を入手することができる。平滑化アクションおよび微分器アクションが別々に図示されているが、平滑化微分器（Ｓａｖｉｔｚｋｙ−Ｇｏｌａｙなど）が選択される場合には、これらを単一のステップ（図１７の１７２）に組み合わせることができる。大きいサポートが利用できるので、正確な微分演算子が、空間領域の関数の近似の代わりに適用され得る。しかし、望まれる場合には、近似を実施することもできる。好ましくは、図１５の方向ベクトルは、画像座標であるが、他の方向を使用することもできる。その後、平滑化関数および微分器関数の周波数領域表現は、特定の画像特性を強調するために、重み付けプロファイルによって操作され得る（たとえば、二乗１５１されることによって）。たとえば、微分１５１および二乗の効果は、低い空間周波数特徴に対して高い空間周波数画像特徴を（ｗｋ）２倍だけ強調することである。空間領域フィルタまたは周波数領域フィルタの任意の組合せが、設計され、図１６の１７１内で使用され得る。Ｎ−Ｄ  Ｓａｖｉｔｚｋｙ−Ｇｏｌａｙフィルタが、図１６の１７１内で使用され得る空間領域フィルタの例である。
【０１３９】
周波数領域表現の説明される計算、図１５の１５１、図１６の１７１、および図１７の１７２は、部分画像自体の使用を必要としない。したがって、これらのステップは、部分画像の計算の前に実行され得る。これらの項の事前計算は、将来の使用のために複数の可能な関数またはフィルタを計算し、記憶することを可能にする。これは、たとえば、複数の可能な平滑化関数を使用すること、または異なる平滑化関数の効率的な試験を可能にする。図１５の平滑化演算および微分演算のこの分離または分断は、大きいサポートを有する平滑化関数またはフィルタ・カーネルの使用という計算的重荷をも軽減する。というのは、これを部分画像ごとに再計算する必要がないからである。上記の方法では、平滑化演算および微分演算は、アーティファクトを克服するために正確な周波数領域一次微分演算子を使用して、別々のステップに分断され、このアーティファクトは、微分器の空間領域表現によって引き起こされる。アーティファクトは、リップルまたは制限された周波数帯応答を含む可能性がある。平滑化微分の周波数領域表現は、その代わりに、平滑化微分カーネル（たとえば、Ｓａｖｉｔｚｋｙ−Ｇｏｌａｙ平滑化微分器）に離散フーリエ変換を適用することによって直接に計算され得る。
【０１４０】
本方法の実施形態は、図１６および図１７の特徴抽出器関数またはフィルタを含むことができる。特徴抽出器関数は、平滑化関数を組み込むか含み、これによって援助され、これの一部であり、またはこれを置換する。フィルタ抽出器関数は、非情報成分（雑音およびドリフトなど）の影響を低減しながら画像情報（画像レジストレーションに有用な）を強調する。最善の関係を達成する方法または最善のフィルタは、画像の空間特性および周波数特性に依存する。いくつかの実施形態では、関係は、部分画像の特性に依存する可能性があり、フィルタは、画像の異なる部分画像について更新され得る。通常の実施形態では、特徴抽出器フィルタは、低周波数のデータおよび高周波数の雑音のＤＣ成分の影響を低減する帯域フィルタである。というのは、ほとんどの画像情報内容が、中間周波数で支配的であるが、雑音およびドリフトが、それぞれ高周波数および低周波数で支配的であるからである。他の実施形態では、特徴抽出器フィルタは、画像に適用される空間領域または周波数領域での帯域フィルタとすることができる。特徴抽出器フィルタの周波数仕様は、画像の有用な特徴を強調し、雑音および画像の重要ではない特徴（ＤＣ成分など）を除去する、Ｎ−Ｄ画像の特性に基づいて定義され得る。適切なフィルタは、望ましくない特徴（雑音およびＤＣ成分など）を除去しまたは改善すると同時に、画像の十分な有用な内容（特徴）を抽出（強調）することのできるフィルタである。このフィルタは、当初に周波数領域の空間領域で定義されるが、通常は、周波数領域の画像に適用される。
【０１４１】
空間領域では、有限差分（一次またはより高次）、窓関数（Ｈａｎｎ、Ｈａｍｍｉｎｇ、Ｂｌａｃｋｍａｎ−Ｈａｒｒｉｓ、ｃｏｎｆｉｎｅｄ  Ｇａｕｓｓｉａｎ、Ｂｌａｃｋｍａｎ−Ｎｕｔｔａｌ、Ｄｏｌｐｈ−Ｃｈｅｂｙｓｈｅｖ）、Ｓａｖｉｔｚｋｙ−Ｇｏｌａｙ平滑化微分器（一次またはより高次）、またはこれらのフィルタの組合せを含む様々な特徴抽出器フィルタが使用され得る。
【０１４２】
周波数領域では、特徴抽出器フィルタは、理論的一次導関数（ｉ×ω）、理論的一次導関数と窓関数の周波数表現との積または組合せ（（ｉ×ω）×（窓関数の周波数表現））、および／または平滑化微分器など、画像情報を強調し、雑音を低減することのできる、カスタム設計の周波数領域関数（ω関数）を含む。
【０１４３】
たとえば、上で議論した平滑化フィルタおよび／またはＳａｖｉｔｓｋｙ−Ｇｏｌａｙ微分器フィルタは、特徴抽出器フィルタである。微分器フィルタは、画像の高周波数情報を強調し、画像の強度勾配を抽出し、平滑化フィルタは、雑音を除去し、より正確なレジストレーションにつながる。しかし、必ずしも微分器フィルタとして定義されるのではない他の帯域フィルタが、同様のプロパティを有する画像に適用され得る。
【０１４４】
図１５の画像１４０に戻ると、部分画像１４１は、図１４を参照して説明したものに似た形で選択される。図１４とは異なって、Ｈａｍｍｉｎｇ窓などの窓関数１４４が、その後、部分画像１４１に直接に適用される（たとえば、平滑化関数なしで）。ＤＦＴ  １４５などの周波数領域変換が、その後、部分画像１４６の間の相互相関を形成するために組み合わされ得る周波数領域部分画像を作るために適用される。多次元画像に関して、Ｎ次元のそれぞれに１つのＮ個の実数ＦＦＴが適用され得る。その後、周波数領域相互相関が、周波数領域平滑化微分器関数表現１５１と組み合わされまたはこれと共に操作されて、周波数領域微分相互相関部分画像を作ることができる。それは、特定の方向での画像の勾配の尺度である。上で説明したように、これは、画像の強度または他の特性に基づくものとすることができる。図１４に関して議論したように、逆変換は、相対変位がそれから計算され得る画像領域相互相関１４７の計算を可能にする。
【０１４５】
ここで図１５、図１６、および図１７に示された完全なプロセスを再検討すると、平滑化関数、Ｎ−Ｄフィルタ、または平滑化微分カーネルは、画像データに直接には依存しないので、事前に計算されまたはライブラリから呼び戻され得る。必要ではないが、平滑化および微分が、別々の演算に分割されることも可能である。おおまかに言って、平滑化は、平滑化関数または平滑化カーネルに周波数変換を適用することによって周波数領域に移動される。この変更のゆえに、平滑化カーネルのサポートのサイズを最小化することに計算的利点はほとんどなく、平滑化関数、Ｎ−Ｄフィルタ、または平滑化微分カーネルのプロファイルの選択におけるより多くの自由と、より少ないリップルを有し（すなわち、より少ない雑音が追加される）、様々な周波数応答に適合され得る、大きいサポートを有するフィルタを使用する可能性とが与えられる。この特定の例では、二乗された周波数領域微分カーネルが、その後、画像のフーリエ変換が計算された後に使用される。部分画像畳み込みパイプラインからの平滑化関数または平滑化カーネルの計算の分断は、計算コストの大幅な節約を実現することができ、計算的に高価な畳み込み演算の必要を除去する、より単純でより高速のアルゴリズムにつながる。さらに、この分断は、画像から有用な画像情報を抽出することによって、他のコンピューティング方法と比較して大幅により高いレジストレーションの精度を達成する方法を提供する。
【０１４６】
図１５、図１６、および図１７の実施形態では、平滑化関数および平滑化フィルタは、画像領域と周波数領域との両方で表現されている。これは、計算の複雑さに対する影響をほとんどまたは全く伴わずに、平滑化演算および微分演算の特性を調整するためのはるかに大きい自由を与える。たとえば、二乗された周波数領域平滑化微分の使用において。周波数領域平滑化微分カーネルの二乗と同等の、微分両方画像の適用の効果は、低空間周波数特徴に対して高空間周波数画像特徴を（ｗｋ）２倍だけ強調することである。我々は、方法が任意の重み付けを適用することを可能にすることに留意する。図１４の方法に関してまたは上の説明で議論した変形は、適用可能な場合に図１５の方法に対しても行われ得る。たとえば、強度の代わりに異なる画像特性が使用され得る。我々は、用語画像が使用されたが、これが、２Ｄ画像と、より高次またはより低次の画像または入力との両方を参照するものと解釈されなければならないことにも留意する。
【０１４７】
図１６の実施形態では、平滑化微分器関数（１５１）は、Ｎ−Ｄフィルタ（１７１）のより一般的なフォーマットで表現される。フィルタは、空間領域または周波数領域で設計され得、単一のω領域Ｎ−Ｄフィルタを形成するために組み合わされ得る。ω領域Ｎ−Ｄフィルタは、画像のプロパティに基づいて選択され得、１からＮまでの任意の数とすることができる。
【０１４８】
図１８は、本方法のワークフローが、少なくとも、副画素部分を任意の個数の次元（２Ｄ、３Ｄ、その他）を有する画像に拡張することを可能にするように適合される、副画素部分での本発明の広義の形態の実施形態を示す。２つの画像が、まず、１画素未満以内でレジスタされ（１７３）、Ｎ−Ｄ窓関数が画像に適用される（１７４）。ＤＣ成分が除去され（１７５）、データが正規化された（１７６）。Ｎ−Ｄ  ＦＦＴが、両方の画像に適用される（１７８）。一方の画像のＦＦＴが共役にされ、位相角が計算された（１７９）。位相値が、エイリアシングおよび／または雑音の影響を打ち消すために境界から除去された。次のステップでは、位相勾配／ｓが、最小二乗法を使用して計算された（１８１）。雑音の多いスプリアス信号を除去するために、±π以上の位相誤差を有する値が、位相データから除外され、位相勾配が、もう一度計算された。このプロセスは、すべての誤差値が±π以内であったか、または反復回数が反復のある個数に達する時まで継続された（１８１）。相対副画素変位が、位相勾配データから計算された（１８３）。
【０１４９】
図１９および図２０は、副画素部分での本発明の広義の形態の第２および第３の代替実施形態である。雑音の多いスプリアス信号は、図１８の１６９に似た手法で除去された。位相データから異常値を除去する余分なステップが、図１９の１８４で使用され、ここで、π／２未満の位相誤差を有する位相値だけが、位相勾配を推定するために選択された。このステップは、π／２を超える位相誤差を有する異常値を除去するのを助ける。異常値は、さらなるステップで、メディアン・フィルタの使用が図２０の１８５および１８６で使用された  ことによって除去された。メディアン・フィルタは、位相データからバースト雑音を除去するのを助け、これは、勾配のより正確な推定値を有するのを助ける。
【０１５０】
図１は、２つの画像の間の並進シフトを見つける本発明の２Ｄバージョンの実施形態を示す。並進シフトは、画像の間の任意の座標系内とすることができる。この流れ図は、副画素精度を有する並進シフトを見つける２ステップ・プロセスを実証するものである。本方法は、低い計算的複雑さを有する２ステップ法を含む。本方法は、多くの既存の方法が上で説明した問題を有する低特徴画像または雑音の多い画像においてさえ、副画素シフトを測定することができる。同一のサイズの２つの画像を入手した１後に（さらなるステップが、一致するサイズを有するようにするために画像をクロッピングしまたは縮小することを含む場合がある）、整数シフトおよび副画素シフトが、別々に見つけられる。整数シフトが、まず計算される２。整数シフトのすべてまたは少なくともほとんどが、１画素精度以内、好ましくは１／２画素精度以内で計算されることを保証するために、頑健な方法が使用される。その後、一方の画像が、シフトを副画素シフトに縮小するためにシフトされる３。好ましくは第１の方法と異なる、第２の方法が、副画素シフトを判定するのに使用される４。整数画素シフトが計算され、その結果、画素シフトが０．５画素未満になる場合には、位相がアンラップされ、全般的に連続した曲線（２次元では２Ｄ曲線）を提示した。この曲線は、好ましくは平面としてモデル化される。方法４は、好ましくは、周波数領域の位相差データを使用する。全体的なシフトは、整数シフトと副画素シフトとを組み合わせる（たとえば、加算する）ことによって見出され得る５。
【０１５１】
本発明の実施形態は、大小両方の合成シフトおよび合成画像回転に関して有効であることが示された。シフトを見つける際の本方法の精度は、１画素の１万分の２〜３程度とすることができる。回転試験では、本方法は、回転された画像内の変位を見つけるための匹敵する技法より性能が優れている。本方法は、画像がガウス雑音または画像システム内の他の一般的な雑音によって劣化した時に、改善された頑健性を提供することもできる。本発明の実施形態では、システムのパラメータは、高精度のレベル、低い複雑さのレベル、および雑音に対する頑健性のレベルの間でトレード・オフするために変更され得る。これは、精度と速度との両方が重要である応用例での、微細な副画素画像レジストレーションまたは変形測定のための使用を可能にする。この方法は、非剛体変換に関する粗から微細までの画像レジストレーション技法に組み込まれることも可能である。本方法は、少数の特徴を有する画像を用いる応用例に特に関連し、あるいは、精度およびほぼリアルタイムの分析が重要である場合には、小から大までの並進シフトを見つけることのできる頑健な副画素レジストレーション方法が要求される。
【０１５２】
剛体変換は、回転および並進からなる。非剛体変換は、異なる画像部分での伸縮と組み合わされた並進および回転を含む。説明される方法は、特に剛体画像レジストレーションのために（特に並進を見つけるために）開発された。しかし、本方法は、非剛体の粗から微細までの画像レジストレーションにも組み込まれ得る。したがって、２つの画像が、まず非剛体レジストレーション方法によってレジスタされ得、その後、我々のアルゴリズムが、副画素レジストレーション部分を実行するのに使用され得る（精度を高めるために）。
【０１５３】
本方法のさらなる利点は、ＧＣが、平均画像強度（０周波数成分は０）に対する依存を除去し、本発明を画像強度の変化に対して相対的に免疫性にすることである。ＣＣまたはＧＣの２Ｄ出力内のピークの近傍に関数をあてはめることによって副画素シフトを見つけることが可能である。これらの方法は、２つの画像の間に非常によい一致を有することに頼り、その結果、関数が、正確にあてはめられ得るようになる。雑音の存在は、このあてはめの精度を低下させる。実際には、強度値変化に起因して、よい一致はめったに発生しない。
【０１５４】
本発明の実施形態では、整数計算ステップ２は、広範囲の画像にまたがって頑健になるように選択される。すなわち、整数ステップ２は、詳細に対するより低い精度を有する可能性があるが、広範囲の画像にまたがって正確な整数シフトを計算することができる。画像は、たとえば、高いまたは低い特徴密度、特徴形状、または雑音レベルを有することによって変化する可能性がある。正確な整数計算は、方法によって使用する  ２つの画像が第１のステップで画素の半分以内にレジスタされることを可能にする。画像が１画素（または１画素の半分またはさらに）以内にレジスタされたので、副画素ステップ４は、情報に関するある種の仮定を行うことができる。具体的には、位相データが副画素ステップ４で使用される場合に、位相データは、アンラップされる必要がなく、直接に処理され得る。既存の方法は、たとえば、特に少数の特徴を有する画像内で苦闘する計算要件（たとえば、相互相関または正規化された相互相関を使用する）のゆえに、頑健な整数シフト法を使用してこなかった。この方法は、整数ステップまたは単一のステップが、非常に正確なレジストレーションを入手することを試みる場合にも低速化する可能性がある。これは、副画素推定値の全体的な精度を低下させるという影響を有する。本発明の実施形態は、少数の特徴を有する画像内であっても整数シフトを見つける際に頑健である整数部分を含む。我々の方法のこの利点は、画像特徴がしばしば豊富ではない多くの実用的応用例にとってクリティカルである。
【０１５５】
頑健性は、雑音の多い画像または低特徴を有する画像内で、セットされた値、好ましくは１画素の値、より好ましくは１画素もしくは０．５画素、またはより小さい画素の値未満の誤差を有する整数シフトを計算する能力と定義することができる。提案される方法での整数シフトを見つけることは、２つの画像特徴の照合に基づく。以前の方法（たとえば、相互相関）は、強度値を照合し、この強度値は、低特徴画像では制限された分散を有し、したがって、よい一致を見つけることは難しい。説明される方法は、ＳＧＤなどの勾配ベースの方法を使用して、まず強度値の勾配を抽出し、その後、ＣＣを使用して勾配値を照合する。勾配法は、好ましくは、平滑化および／または形状保存特徴を有する。形状保存特徴は、好ましくは、勾配の基礎になる真の値を最小限に変更する。したがって、説明される方法は、低特徴画像からより多くの特徴を抽出することができ、これが、照合プロセスをより正確なものにする。
【０１５６】
本発明のさらなる実施形態では、本方法は、
異なるウィンドウイング関数と、
位相差データの正規化および事前フィルタリングと、
エイリアシングの影響を減少させ、シフト推定値の精度を向上させるための変更と
のうちの任意の１つまたは複数と共に画素シフト計算（好ましくは、位相ベースの方法を使用する副画素シフト計算）を組み込むことができる。
【０１５７】
図２は、整数シフト・ステップ２の実施形態を示す。特定の実施形態では、Ｓａｖｉｔｚｋｙ−Ｇｏｌａｙ微分器が、勾配相関を入手し、整数シフトを見つけるのに使用される。本発明の実施形態は、２つの画像の間の並進シフトを見つけるために、２Ｄ−ＣＣ（２Ｄ相互相関）への強度勾配入力を使用する。本方法が、ここでは強度勾配に関して説明されるが、他の画像特性が、実際面で使用され得、あるいは、特性が、画像レジストレーションが行われる前に変更され得る。一般に、特性値は、画像の画素のそれぞれから入手されるが、ある種の場合には、異なるまたは特定の複数の点が選択される。Ｉ１（ｘ，ｙ）およびＩ２（ｘ，ｙ）が、点（ｘ，ｙ）でのＮ×Ｎのサイズを有する２つのグレイスケール画像の強度値を表すものとすると、２つのテンプレートの２Ｄ−ＣＣは、
【数２】
2019517073.tif 000004
と定義され、ここで、−Ｎ＜ｋ、ｊ＜Ｎであり、上線は、複素共役を表す。負のインデックスは、通常、０に置換される（すなわち、０パディング）。いくつかの実施形態では、２Ｄ  ＦＦＴおよびその逆（２Ｄ  ＩＦＦＴ）を使用して周波数領域で２Ｄ−ＣＣを計算することが、より簡単である。
【数３】
2019517073.tif 000005
Ｉ１およびＩ２は、画像の強度行列であり、Ｎ×Ｎのサイズを有する画像の２Ｄ  ＦＦＴは、
【数４】
2019517073.tif 000006
によって与えられ、ここで、（ｉ＝√−１）である。いくつかの実施形態では、正規化されたクロスパワー・スペクトルが、並進シフトを見つけるのに使用される。
【数５】
2019517073.tif 000007
【０１５８】
この実施形態は、並進シフトを見つけるために画像強度を強度勾配（一般に複素項の形で記述される）に置換することによって勾配相関（ＧＣ）に適合され得る。好ましい実施形態では、実数部は、第１の方向での強度勾配に関し、複素部は、第２の直交方向での強度勾配に関する。この複素項の実数部および虚数部を計算するために、それぞれ水平および垂直の強度勾配が近似される。代替技法は、勾配を計算するのにｘ方向（ＣＤｘ）およびｙ方向（ＣＤｙ）での中心差分を使用する。
ＣＤｘ＝Ｉ（ｘ＋１，ｙ）−Ｉ（ｘ−１，ｙ）
ＣＤｙ＝Ｉ（ｘ，ｙ＋１）−Ｉ（ｘ，ｙ−１）
【０１５９】
中心差分は、特定の画素位置での理想的な微分を近似する。一次中心差分および二次中心差分が使用されてきた。好ましい実施形態では、本方法は、多項式への移動最小二乗あてはめを使用して信号を平滑化することによって理想的な微分器のより正確な近似を全般的に提供するＳａｖｉｔｚｋｙ−Ｇｏｌａｙ平滑化微分器（ＳＧＤ）を使用する。これは、特に雑音の多い信号に関して、中心差分より頑健であり、より正確である。というのは、ＳＧＤが、雑音による近傍点の破損によってそれほど簡単には影響されないからである。すなわち、補間は、近傍点に直接には依存せず、平滑化効果を有する。補間が近傍点に直接には依存しないので、精度は、ＳＧＤの精度を高めることによって高められ得る。
【０１６０】
ＳＧＤの使用に起因するさらなる改善は、信号形状の保存である。信号形状の保存または維持は、信号から雑音を除去することと、近傍点から情報を収集することとを含むことができる。すなわち、信号の全体的な形状ならびに雑音低減が重要である。これは、推定値を作成するために信号点の両側の複数の点を使用することによって達成され得る。ＳＧＤのさらなる有用な特徴は、適当なフィルタ長の選択によって、信号導関数の分解能を維持できることである。理想的なデジタル微分器は、高周波数雑音を増幅するという不利益を有する。ＳＧＤは、低周波数では理想的な微分器の周波数応答に近く、信号形状を保存するが、より高い周波数では雑音の影響を減衰させる。したがって、一実施形態では、ＳＧＤは、信号形状を保存でき、雑音を除去しながら実際の信号データを保存しまたは最小限に変更することを試みる低域通過微分器とみなされ得る。これは、高周波数データの除去を含むことができる。というのは、雑音が、しばしば、信号の素早い変化として現れるからである。したがって、形状保存機能または形状保存効果は、信号を平滑化する。除去される雑音の量にトレード・オフがある。というのは、これがデータの消失を引き起こすからである。ＳＧＤなどの関数は、画像強度勾配を抽出する微分器を含み、真の画像データを保存しながら雑音を除去することもできる。
【０１６１】
特定の実施形態が、勾配情報を抽出するためのＳＧＤを使用して説明されたが、この手法は、より幅広い様々な方法を使用して実施され得る。代替の方法によって達成され得る、ＳＧＤ勾配推定法の１つの利点は、中心差分ベースの方法と比較して異なる特性をもたらす、より高い周波数およびより低い周波数で異なる周波数応答を有することである。図４は、中心差分およびＳＧＤの周波数応答４０、４１、および４２を示し、両側に周波数帯を有する、約０．７の正規化された周波数に０減衰または高い減衰がある。この周波数応答は、雑音を除去し、信号形状を保存し、低周波数で理想的な微分器に匹敵することができる。すなわち、この周波数応答は、整数画素レジストレーションに有用なデータである信号対雑音比を向上させることができる。周波数応答は、好ましくは、特に低周波数で、信号形状または真の基礎になる値を保ち、特に高周波数で、信号から雑音を除去するために帯域フィルタまたはノッチ・フィルタに似る。説明される特定のシステムに対する可能な代替案は、Ｌａｎｃｚｏｓ微分器または異なる次数のＳＧＤの使用を含む。
【０１６２】
本発明の一実施形態では、勾配推定値は、離散信号に関して、それらが、移動最小二乗あてはめを使用するのではなく、好ましくはフィルタ次数ごとの係数の表を用いて、畳み込み法を使用して簡単に実施され得るという利点を有する。これは、計算の複雑さを低減することによって本方法の計算コストを低下させ、ハードウェア実施（たとえば、プロセッサ、ＦＰＧＡ、またはＧＰＵ内での）を実現可能にする。畳み込み法は、複数の勾配推定方法に関して、また、複数の次数のエスティメータに関して使用可能である。特定の実施態様では、７点三乗一次導関数ＳＧＤが使用される。これは、周波数応答プロパティに基づき、近傍点からのより多くの情報を使用する。このフィルタの畳み込みカーネルは、
ＳＧＤ（カーネル）＝［２２，−６７，−５８，０，５８，６７，−２２］）／２５２
である。
【０１６３】
図４は、このカーネル４０を示し、このカーネルの周波数応答を一次中心差分４１および二次中心差分４２と比較する。一次中心差分および二次中心差分の周波数応答は、周波数帯全体にわたって類似するが、ＳＧＤは、正規化された周波数の半分より高い周波数４３に関して異なる応答を有する。
【０１６４】
特定の実施形態では、以下の実施プロセスに従うが、当業者は、本明細書で説明する項またはステップの組合せにおいて変形形態が可能であることを理解するであろう。ＳＧＤの畳み込みカーネルは、それぞれＳＧＤｘおよびＳＧＤｙを形成するために、画像の各行および各列に適用される２０。これらは、我々の方法の複素項を形成するために組み合わされる。
ＳＧＤ（ｘ，ｙ）＝ＳＧＤｘ（ｘ，ｙ）＋ＳＧＤｙ（ｘ，ｙ）×ｉ
【０１６５】
この項は、２Ｄ  Ｓａｖｉｔｚｋｙ−Ｇｏｌａｙ勾配相関（ＷＧＧＣ）
【数６】
2019517073.tif 000008
を形成するのに使用され２１、ここで、ＳＧＤ１およびＳＧＤ２は、長方形行列の形になっている。
【０１６６】
図５ａは、元のパワー・スペクトル５０を示し、図５ｂは、雑音が存在する、同一のスペクトル５０を示す。図５は、通常、イメージング・システムの制限された解像度に起因して、素早い強度変化がほとんどの画像内でまれであることを実証するものである。０周波数５２（すなわち、ＤＣ成分）は、画像の中心へシフトされ、より高い周波数は、画像の周辺部分５３または外側により近い。画像の微細な詳細は、通常、雑音成分の周波数と比較して、高周波数成分のより低い範囲内にある。白色ガウス雑音は、図５（ａ）の画像の微細な詳細と比較して、図５（ｂ）内でより高い周波数に広がっている。これは、図５ｂの周辺付近の精細度の消失に見られ得る。
【０１６７】
図６は、複数の窓関数を示す。本発明の実施形態では、窓関数は、勾配画像に適用される２２。本発明のさらなる態様は、周波数領域技法を使用し、複数の窓関数とパラメータおよび／または調整とが考慮される。窓関数６０のパラメータは、好ましくは、それぞれ整数部分または整数ステップ２および副画素部分または副画素ステップ４に適当になるように選択される。窓関数６０は、境界効果、エイリアシング、および雑音を減少させることができる。これは、窓関数が、長方形の窓によって引き起こされるスプリアス高周波数成分の導入を防ぐからである。雑音と画像の微細な詳細（エッジなど）との両方が、周波数領域の高周波数帯にある。したがって、適当な窓関数の選択は、同時に微細な詳細を保存し（２つの画像の間の正確な照合を可能にするため）、雑音を低減するために重要である。
【０１６８】
通常、窓関数は、画像情報の保持と雑音、エイリアシング、および境界効果の除去との間にトレード・オフを有する。したがって、窓関数（ＢｌａｃｋｍａｎまたはＴｕｋｅｙなど）は、以前には、雑音を除去するなんらかの方法の中で使用された。しかし、これらの窓関数は、微細な画像詳細をも除去し、これが、より不正確なシフト推定を漏らす。アルゴリズムが雑音に対して頑健である（たとえば、形状保存勾配推定の使用に起因して）場合には、Ｈａｍｍｉｎｇ窓など、より低減衰の窓が使用され得る。これの利点は、画像特徴がより詳細に維持されることである。すなわち、頑健な整数勾配エスティメータの使用は、信号対雑音比を改善するために改善された窓関数と組み合わせられ得る。
【０１６９】
異なる窓関数は、異なる特性を有する。いくつかの既知の窓関数は、画像から雑音を除去するのを助けるが、高周波数成分（画像の微細な詳細）をも除去する。既存の方法は、雑音に対処するために２つの窓関数を使用してきた。しかし、窓関数のこの選択は、画像の微細な詳細を除去するので、精度を低下させる。その代わりに、高周波数成分を保存する窓関数が選択され、他の技法が、雑音に対処するのに使用される。整数レジストレーションまたは副画素レジストレーションに関して提案される可能な技法は、正規化アルゴリズムの改善、メディアン・フィルタの使用、および２Ｄ回帰の前に位相データを選択するための関数の使用を含む。いくつかの場合に、実質的に同一のまたは類似する技法が、副画素シフトと整数シフトとの両方に使用され得る。
【０１７０】
特定の実施形態では、周波数応答が減衰の穏やかな傾き６０を有し（−６ｄＢ／オクターブ）、第２ローブ６２内で高い減衰（−４３ｄＢ）を有するので、Ｈａｍｍｉｎｇ窓が使用される。これは、整数シフトに関して、低周波数データの大きい部分がよく保存される（穏やかなロールオフに起因して）と同時に、すべての雑音（高周波数に集中する）がよく減衰されることを意味する。その後、窓関数は、勾配推定法の雑音頑健性プロパティに関して選択される。たとえば、Ｓａｖｉｔｓｋｙ−Ｇｏｌａｙ勾配相関（ＳＧＧＣ）は、雑音の存在下であっても、画像の微細な詳細を抽出するのを助ける。これは、ＳＧＤカーネルまたはＳＧＤ法が、相関が計算される前の信号から信号雑音を除去し、低減し、または改善するからである。改善された雑音除去技法は、より多くの信号（真の画像）データが保たれることをも可能にする。勾配推定法と窓関数とのこの組合せは、雑音の低減とより高い周波数での画像の詳細の維持との間の妥協を提供する。比較して、Ｂｌａｃｋｍａｎ窓は、第２ローブの後に激しく減衰し、画像内の高周波数情報（すなわち、微細な詳細）のほとんどを除去し、Ｔｕｋｅｙ窓は、その第２ローブ内で十分な減衰を有しておらず、多すぎる雑音が残ることをもたらす。Ｔｕｋｅｙ窓は、その周波数応答のリップルに関して位相ひずみをも生じる。
【０１７１】
図７は、改善された正規化方法２３を示す流れ図である。この正規化手順は、本方法のダイナミック・レンジおよび雑音に対する頑健性を改善することを目標とする。上で説明した正規化されたクロスパワー・スペクトル（ＮＣＰＳ）およびＮＧＣの個々の項は、ＤＣ成分を含んだ。ＤＣ成分の包含は、これが信号の大きさと周波数のそれぞれの相対重みとに影響したので、信号のダイナミック・レンジの低減をもたらす。ＤＣ成分は、平均画素値を反映する。図７の流れ図７０は、ＳＧＤ１およびＳＧＤ２が別々に正規化される（最大絶対値７２（または推定された最大値）によって除算する平均（ａｖｅｒａｇｅ）（たとえば平均（ｍｅａｎ））値７１を減算することによって実施形態を示す。最大値は、通常は平均値の除去の後に計算される。その後、正規化された値７３が、ＳＧＧＣ式に適用される。これは、ダイナミック・レンジを増大させ、強度感度を減少させるという利点を有する。これは、複素項の実数部および虚数部が、平均値の減算によるＤＣ除去の後に別々に正規化され、お互いまたは画像平均値によって影響されないからである。代替実施形態では、ＤＣ成分は、たとえば０の周波数成分の振幅を入手することによって、周波数領域で計算され得る。
【０１７２】
図３は、副画素整数シフト・ステップ４の実施形態を示す。特定の実施形態では、変更された位相ベースの方法が、副画素シフトを見つけるのに使用される。図１を参照すると、好ましくは、副画素レジストレーション４は、画像が１画素未満以内、好ましくは０．５画素未満以内でレジスタされた後に実行されるが、いくつかの状況で（たとえば、整数シフトが既知である場合または存在しないことが既知である場合に別々に使用され得る。副画素シフトは、好ましくは、周波数シフト手法に基づいて計算される。すなわち、この計算は、周波数領域で、位相差などの位相特性の勾配を識別することによって副画素シフトが計算され得るという特徴を使用する。前と同様に、システムは、通常、強度値を用いて評価されるが、要求される場合に、ある範囲の画像特性が使用可能である。使用される複数の点は、好ましくは、整数シフトに関して使用されるものと同一であるが、異なる複数の点が使用され得る。
【０１７３】
信号Ｆ（ｗ）のフーリエ変換は、別々の実数部および虚数部を含み、次の形で表現され得る。
【数７】
2019517073.tif 000009
２つの信号の位相差（φ）は、次式によって与えられる。
【数８】
2019517073.tif 000010
【０１７４】
本発明の実施形態では、２次元データ信号の位相差（φ）データの傾きが、副画素シフトを見つけるのに使用される。これは、位相ラッピングが行われなかった時、たとえば、すべての整数シフトが計算され、除去された場合に、相対的に単純である。位相差の傾きを計算するステップは、画像内の副画素並進シフトを見つけるために、２Ｄおよび位相面（φｐ）に拡張可能である。
【０１７５】
本発明の実施形態では、本方法の整数部分に関して前に説明した正規化手順７０が、副画素画像に適用され得る３１。好ましい実施形態では、窓関数が使用される３０。好ましい窓関数は、図６に示されたＨａｎｎ窓関数である。これは、位相差を計算する３３に使用される強度値のフーリエ変換（好ましくはＦＦＴを用いる）を計算する３２前に画像に適用される３０。Ｈａｎｎ窓は、周波数領域でのそのプロパティのゆえに選択された。Ｈａｎｎ窓周波数応答は、第２ローブ内での３１．５ｄＢの適度な減衰および減衰の急な傾き（−１８ｄＢ／オクターブ）を有する。この減衰は、たとえばＨａｍｍｉｎｇ窓と比較して、第２ローブ内ではより小さいが、より急な減衰の傾きを有する。副画素ステップでは、画像は、既に１／２画素以内でレジスタされている（より大きい差は、整数部分でレジスタされ、シフトによって除去される３）。したがって、副画素ステップは、より大きい特徴ではなく雑音の除去に焦点を合わされている。これは、Ｂｌａｃｋｍａｎ窓関数およびＨａｎｎ窓関数が雑音の除去に使用されなければならないことを示唆する。というのは、これらの両方が、高周波数データ内の雑音をフィルタリングできるからである。しかし、Ｈａｎｎ窓は、Ｂｌａｃｋｍａｎ窓と比較して、その第２ローブ内でより少ない減衰を有するので、相対的に低い周波数のデータをも保存する。すなわち、第１ローブおよび第２ローブの相対減衰は、Ｂｌａｃｋｍａｎ窓の場合がより大きい。
【０１７６】
図１１は、スプリアス周波数がデータに導入される場合の、信号１１１、図１１（ａ）に対するエイリアシングおよび／または雑音の可能な影響１１０、図１１（ｂ）を示す。位相データからの副画素シフトの推定を改善するために、一般に雑音またはエイリアシングによって引き起こされるこれらのスプリアス周波数が、除去され得る。本方法の実施形態は、スプリアス周波数を除去するための、データに対する２Ｄメディアン・フィルタを含む。メディアン・フィルタのパラメータまたは使用されるフィルタリング方法が変化する可能性があることを理解されたい。他のデジタル・フィルタリング方法（たとえば、平滑化フィルタ、平均化フィルタ、およびガウシアン・フィルタ）およびフィルタ半径が、適切な実施態様を見つけるために変更され得る。特定の実施形態では、２Ｄメディアン・フィルタの近傍サイズが、相対的に小さくなるように（すなわち、３×３）選択された。これは、中心での傾きを変更する（たとえば、測定される傾きに対する影響を低減する）ことなくデータを平滑化するのを助ける。この２Ｄメディアン・フィルタは、傾きまたは勾配を見つけ、副画素シフトを計算するのに使用される２Ｄ回帰（または他の方法）に使用されるデータを改善するという効果を有する。２Ｄ回帰は、データが雑音によって影響される時に、より不正確になる。
【０１７７】
図１３は、メディアン・フィルタに対する雑音の影響を低減する方法の流れ図である。これは、近傍内のデータの中央値を選択することと、推定された平面をデータに適用することとによって、変動性を低減しまたは位相データを平滑化する。関連する影響は、代替の平滑化手段によって達成され得る。他の実施形態では、メディアン・フィルタは、より高い周波数で最も顕著な雑音を目標とする非線形フィルタによって置換され得る。複数の、たとえば第１および第２の、画像が、まず選択され３１、位相データ３２が、シフトを解釈するために計算され得る。その後、メディアン・フィルタ３３が、位相データを平滑化するのに使用され得る。次のステップでは、周波数制限または範囲の縮小が、信号に適用される３４。雑音の量を判定する３５ことと、その後、適当なサンプル・サイズを計算しまたは入手する３６（例えば、雑音の量に基づいて）こととによって、より明瞭な位相差が、望まれない周波数を除去する３７ことによって使用され得る。通常、これは、雑音周波数が一般にスプリアス高周波数なので、帯域タイプ・フィルタまたは高域タイプ・フィルタとして実施される。２Ｄ回帰または他の勾配計算技法が、位相差の傾きを計算し、これから副画素シフトを識別するのに使用され得る３８。
【０１７８】
以前の方法は、サンプリングされるデータに関する固定値（たとえば比率）を選択することによって、原点付近のデータの周波数範囲を制限してきた３４。たとえば、２Ｄ回帰で使用されるサンプルの個数（Ｎφ）としてデータの中心付近の画像サイズの０．６という固定値を選択することであり、ここで、値は、試行錯誤によって選択される。本発明の実施形態では、サンプルの適当な個数（Ｎφ）は、画像の雑音レベルの特性に依存する３５。画像の微細な詳細と雑音との両方が、周波数領域データの高周波数成分に存在するので、サンプルの個数（Ｎφ）を正確に選択することが重要である。すなわち、データ中心付近のより小さいＮφを選択することによって、雑音および微細な詳細が同時にデータからフィルタリングされ、微細な詳細の低減およびより低い精度を引き起こす。
【０１７９】
改善された方法の特定の実施形態  データの中心付近のサンプルの異なる個数Ｎφが選択される３７。たとえば、画像が低いレベルの雑音を有する場合には、多数のサンプルが、画像の微細な詳細を最もよく保存するために使用され得る。その代わりに、高いレベルの雑音が存在する場合に、サンプルの個数すなわちサンプル・サイズＮφは、このバランスをとるために低減され得る。画像の雑音レベルまたは画像内に存在する雑音の特性が、画像雑音推定の当業者に既知の複数の方法によって計算されまたは判定され得る３５。たとえば、データの２Ｄ標準偏差（２ＤＳＤ）が使用され得る。統計モデル化を使用して信号雑音を推定する、より洗練された形がある。これらは、雑音推定技法またはブラインド雑音推定（ｂｌｉｎｄ  ｎｏｉｓｅ  ｅｓｔｉｍａｔｉｏｎ）を含む。２ＤＳＤは、データ内の雑音のバランスのとれた見方を提供するので実用的である。雑音およびエイリアシングがないデータという理想的な場合には、２ＤＳＤは、完全に対称であり、０の２ＤＳＤ値をもたらす。
【０１８０】
雑音計算法は、データ内に存在する雑音の量を表現する値または特性を生成するのに使用され得る。その後、関数または関係が、雑音レベルに依存して複数のサンプルＮφを定義しまたは入手するのに使用され得る３６。たとえば、線形関係が適用され得る（上限と下限との間で）。この関係の例の形は、次の通りである。
【数９】
2019517073.tif 000011
Ｎφ＝ｐ×画像サイズ        （１２）
【０１８１】
この場合に、Ｎφの最大値は、画像サイズの０．９５であり（境界効果を防ぐために）、最小値は、画像サイズの０．７５である（画像詳細を保存するために）。変数ｐは、画像サイズ（または採用されるサンプル）の変化を補正するのに使用される。最大値および最小値は、雑音、データ・セット、もしくは画像（たとえば、特徴周波数またはサイズ）、または他のパラメータに関してより多くの情報が既知である場合に、変更され得る。他の実施形態では、関係は、複数のサンプルをよりよく計算するために、非線形または複素とすることができる。いくつかの実施形態では、主成分分析（ＰＣＡ）および特異値分解（ＳＶＤ）が、２Ｄ回帰の代わりに、位相面データの勾配を見つけるのに使用され得るが、これらは、通常はより多くの計算を必要とする。
【０１８２】
図８は、合成シフトを適用することによって本発明がそれに対して評価され得る６つの標準ランドサット画像を示し、その後、各方法から生じるシフト値を推定した。たとえば１画素より小さい値または１画素より大きい値からなる２つのグループ内の、様々な副画素シフトまたは整数シフトが、ＦＦＴシフトを使用して適用され得る（これが実世界の実験に似ているので）。図９は、画像８０が、必要な場合にあるサイズまたは部分画像９０にどのようにクロッピングされ得るのかを示す。部分画像９０の中心点は、２０画素のステップ増分と周辺への６４画素の最小距離とを用いて画像全体にまたがって選択された。これは、画像ごとに４００個の部分画像を、６つのランドサット画像８０について合計２４００個の部分画像をもたらした。画像を部分画像に分割することは、狭い地域の、地域の、一部の、またはすべての変位または変形の計算を可能にする。画像全体の変形マップは、部分画像のそれぞれからのデータの組合せである。部分画像の使用は、大きい面積にわたる回転または他の不均一並進変化に関して有利である可能性がある。選択される部分画像のサイズは、変形／変化のタイプおよび画像の固有の特徴に依存して変化する可能性がある。
【０１８３】
並進シフトを見つける際の平均レジストレーション誤差（ＲＭＳＥＴ）は、
【数１０】
2019517073.tif 000012
に基づいて計算されたが、ここで、［ｘ，ｙ］は、画素位置であり、ハット記号は、推定された画素位置値を示し、ｍは、ランドサット画像８０番号であり、ｎは、部分画像９０番号である。レジストレーション誤差または他の誤差は、平均平方の意味で計算され得る。ＲＥＴ値は、並進シフトを見つける際の方法の性能を評価するのに使用された。適当な窓関数が、画像レジストレーション方法の性能を改善するのに使用される。具体的には、既存の方法は、Ｈａｎｎ窓の追加によって改善され得る。
【０１８４】
回転の存在下で並進シフト誤差を評価するために、０．５°、１°、２°、および３°の回転値ならびにが、画像中央の回転中心を用いて６つのランドサット画像のすべての画像に適用された。画像８０は、回転を適用する前および回転を適用した後に、部分画像にクロッピングされた。その後、ｘ方向およびｙ方向の並進シフトが、説明される方法によって計算された。誤差の推定値は、２次元の並進に基づいて計算され得る。これは、部分画像９０ごとに並進シフトの計算された大きさＴｘ、Ｔｙを比較することを含む。
【数１１】
2019517073.tif 000013
【０１８５】
この並進シフトは、画像全体に連続的な回転パターンを与えるために補間され得る。部分画像９０ごとの理想的な回転からのシフトの大きさは、
【数１２】
2019517073.tif 000014
としてモデル化され得、ここで、シータは、ラジアン単位の回転角であり、［ｘ０，ｙ０］は、画像の中心であり、［ｘｒ，ｙｒ］は、回転された画像内の点［ｘ，ｙ］の画素位置である。比較は、回転に関する平均レジストレーション誤差とすることができ、
【数１３】
2019517073.tif 000015
であり、ここで、ｍおよびｎは、画像および部分画像の参照番号である。
【０１８６】
本発明の実施形態は、画像の合成シフトに基づいて試験された。この技法は、画像を試験し、方法のシフト推定誤差を見つけるために、既知の合成シフトを適用する。本方法は、他の当事者によって説明された既存の方法と比較された。白色ガウス雑音が、通常、シフト推定の前に画像に追加される。２つの誤差メトリックが、我々の方法の整数部分および副画素部分の精度を評価するために定義された。これらのメトリックは、特定の画像に適用された時と、剛体並進試験を実行することが不可能である場合との、方法精度の推定値を提供する。
【０１８７】
誤差メトリックまたは精度を識別する他の方法が、整数シフト識別および副画素シフト識別の精度を分析するのに使用され得る。精度の推定は、測定値の信頼メトリックを有することが重要である応用例で特に関心を持たれる。これらのメトリックは、異なる特徴および雑音レベルを有する別個の異なる画像での方法精度を推定する。２Ｄ−ＣＣまたはＧＣの出力内の支配的ピークの弁別性は、しばしば、整数シフトを判定する際の方法の能力を示すと考えられる。整数部分に関して、内の正規化されたＳＧＧＣ値が０．８５より大きかった点の個数が、我々の提案される方法で整数シフトを見つける際の誤差を示すのに使用された。このメトリック値が１に近ければ近いほど、整数部分で本方法によって達成されるシフト推定誤差が小さくなる。
【０１８８】
図１４は、上で説明され、図２に示された実施形態に類似する実施形態を示す。２つの部分画像１４０が、より大きい画像１４１から選択されて示されている。代替の場合には画像全体が使用され得るが、その代わりに、主画像内から選択された複数の部分画像を使用することが有益である場合がある。部分画像の最良のサイズは、精度と局所性との間の相対プレーオフならびに画像の雑音特性によって決定される。画像の情報または特徴サイズは、サイズとオーバーラップとの両方に影響する。正確なオーバーラップを計算するために、部分画像は、共通の特徴の実質的なオーバーラップを有しなければならない。というのは、これが、一致を作成するのに十分な情報を提供するからである。
【０１８９】
その後、部分画像は、Ｓａｖｉｔｚｋｙ−Ｇｏｌａｙ導関数などの平滑化関数１４２によって処理される１４３。平滑化関数１４２は、画像（より高い次元での形状または入力とされ得る）の次元ごとに１回適用され得る。実数フィルタが使用されている場合には、２次元推定値が、単一の複素推定値を形成するために組み合わされ得る。より高い次元では、単一の複素ＦＦＴの代わりに、複数の実数ＦＦＴが、次元ごとまたは画像特徴ごとに使用され得る。ピークが、実数ＦＦＴの大きさに基づいて見つけられ得る。その後、前に説明したように、Ｈａｍｍｉｎｇ窓などの窓関数１４４が、スペクトル漏れを制限するのに使用され得る。Ｈａｍｍｉｎｇが使用されるが、ｃｏｎｆｉｎｅｄ  Ｇａｕｓｓｉａｎ窓、Ｂｌａｃｋｍａｎ−Ｎｕｔｔａｌｌ窓、またはＤｏｌｐｈ−Ｃｈｅｂｙｓｈｅｖ窓を含む他の変形形態が、その代わりに使用され得る。その後、変更された部分画像は、画像領域から周波数領域表現に変換される。これは、好ましくは、離散フーリエ変換（ＤＦＴ）または高速フーリエ変換（ＦＦＴ）などのフーリエ変換１４６を介する。周波数領域変換が説明されるが、当業者は、離散ウェーブレット変換（ＤＷＴ）などの他の領域の変換が存在し、同様の結果を提供できることに気付くはずである。その後、周波数領域相互相関１４７が、変更された第２の部分画像の複素共役との変更された第１の部分画像の乗算によって実行される。計算された後に、逆変換が、相互相関を画像領域に戻って変換し、これから、相対変位が計算され得る。たとえば、相互相関は、画像領域相互相関の要素の最大の複素数絶対値を見つけることによって計算されまたは入手される１４８。
【０１９０】
副画素部分の誤差メトリックは、φｐにあてはめられた平面の平面あてはめ誤差と定義された。たとえば、あてはめ誤差は、ノルムまたは最小二乗法によって計算され得る。小さいあてはめ誤差は、副画素シフトの推定における方法のよい精度を示す。小さい近傍を用いる２Ｄメディアン・フィルタを使用する位相データの事前フィルタリングは、位相データから異常値を除去するのを助け、頑健な平面あてはめをもたらし、したがって信頼できる副画素誤差メトリックをもたらす。
【０１９１】
テーブル２およびテーブル３は、既存の方法の選択物のそれぞれの、ｘ方向およびｙ方向のシフト（［ｘ，ｙ］）を示す。我々の方法の平均値は、１画素未満のシフト（テーブル２）に関して、次によい方法（すなわち、Ｓｔｏｎｅ他、ＵＳ６６２８８４５）より約４倍小さい。これは、我々が副画素変位測定に関して同様の手順を選択したが、我々の変更が精度を大幅に改善したことを実証する。
【０１９２】
テーブル２
副画素シフトに関する２４００個すべてのランドサット部分画像（１２８画素×１２８画素）の平均レジストレーション誤差（画素単位）（ＲＥＴ）が、第１列に示される。最後の行は、我々の提案した方法の平均レジストレーション誤差に関して定量化された平均レジストレーション誤差の相対差を提供する。
【表１】
2019517073.tif 000016
【０１９３】
テーブル３
１画素を超えるシフトに関する２４００個すべてのランドサット部分画像（１２８画素×１２８画素）の平均レジストレーション誤差（画素単位）（ＲＥＴ）。最後の行は、我々の提案した方法の平均レジストレーション誤差に関して定量化された平均レジストレーション誤差の相対差を提供する。
【表２】
2019517073.tif 000017
【０１９４】
テーブル３は、我々の方法の、整数シフトおよび副画素シフトを見つける能力を定量化するものである。Ｖａｎｄｅｗａｌｌｅの方法は、位相ラッピングのゆえに、１画素を超える副画素シフトを見つけることができない。Ｓｔｏｎｅ他の方法（ＵＳ６６２８８４５）を使用した結果は、貧弱なテクスチャを有する部分画像内の整数シフトを見つけることに関してＣＣが有する問題点を示す。いくつかの場合に正しい整数シフトを見つけることに失敗したＧｕｉｚａｒの方法を使用する時に、同様の問題が生じた（度合はより低いが）。Ｆｏｒｏｏｓｈ他の方法は、ＣＣに対する位相補正（ＰＣ）の利点のゆえに、穏当に信頼できるものであった。しかし、これらの結果は、オーバーラップしない領域（領域が一方の画像には存在するが他方の画像には存在しない場合）に関する誤差の増加という問題を強調する。以前に公表された方法のうちで、Ｆｏｒｏｏｓｈの方法は、適切な（この場合にはＨａｎｎ）窓の追加によって改善されることが示された。Ｆｏｒｏｏｓｈは、窓関数の悪い選択が、Ｆｏｒｏｏｓｈ＋ＨＷが最もよく実行した画像情報を失うことをもたらす可能性があるので、窓を使用しないことを選択した可能性がある。本発明方法の実施形態は、最も近い既存の方法より５８倍小さい平均値を提供し、小さいシフトおよび大きいシフト（それぞれ、テーブル２およびテーブル３）の計算に関して一貫して良好に動作した。テーブル２およびテーブル３のＦｏｒｏｏｓｈとＦｏｒｏｏｓｈ＋ＨＷとの間の精度の比較は、精度を高める際の適当な窓関数の役割を強調するものである。
【０１９５】
ガウス雑音は、様々なステージでデータ収集方法またはデータ処理方法に入るランダム雑音に起因して、ほとんどのイメージング・システムでの支配的なタイプの雑音である。我々の方法の雑音に対する頑健性を調査するために、ガウス白色雑音（正規化された強度値に対する０の平均値、分散）が、パリのランドサット画像に追加された。レジストレーション誤差が、４００個の部分画像９０にわたって［３．２５０，４．７５０］の［ｘ，ｙ］シフトの推定に関して計算され、平均をとられた。
【０１９６】
図１０は、本方法が、適用されたガウス雑音範囲のそれぞれに関して既存のシステムと比較してかなりより低い推定誤差を有することを示す。本発明の実施形態の平均レジストレーション誤差は、大量のガウス雑音がある（図１０ａ、雑音は正規化された強度値の分散０．１２を有する）場合であっても、０．７８画素であった。本方法の精度を示す際の整数誤差メトリックおよび副画素誤差メトリックの能力を評価するために、この２つの誤差メトリックが、我々の方法に関して各ガウス雑音レベルで計算された（図１０（ｂ）〜図１０（ｃ））。これらの図は、同一レベルのガウス雑音での、本方法の平均レジストレーション誤差と整数誤差メトリック値との間の直接の関係を実証するものである。図１０ａのレジストレーション誤差の増加の緩やかな傾きは、１から約２までの整数誤差の徐々の増加をもたらした。さらに、図１０ｂは、かなりのガウス雑音の存在下であっても、ＳＧＤまたは類似する技法の固有の雑音頑健性を利用することによって、本発明の整数部分は２つの部分画像９０の間の正しい整数シフトを見つける際に良好に動作することを示す。０．８５を超える正規化されたＳＧＧＣ値を有する点の平均個数は、０．１２の分散（σ２＝０．１２０）を有するかなりのガウス雑音を伴う場合であっても、すべての部分画像９０に関して２点未満であった。
【０１９７】
図１０ｃは、副画素誤差メトリックに関して、値が始めから急速に増加し、本方法の副画素精度が０．０００１画素から０．２画素に減少したことを示す。その後、副画素誤差メトリック値の変化は、副画素精度の減少の速度も低下したので、より小さくなった。副画素誤差メトリック値のこの傾向は、このメトリックが、副画素シフトの推定での本方法の精度の表示になるのに十分に敏感であることを示した。
【０１９８】
微細なレジストレーションのための既存の方法では、レジストレーション方法は、擬似極座標ＦＦＴを使用して回転を推定することができる。しかし、回転推定値は、並進シフト推定値ほど正確ではなく、通常は、回転中心が画像の中心にある時に回転を見つけるのに使用される。ほとんどの実用的な応用例では、回転は、非剛体画像変換の一部である。このタイプの回転は、剛体画像レジストレーション方法の精度を低下させ、擬似極座標ＦＦＴは、それを正しく推定することができない。
【０１９９】
図１２は、ミシシッピのランドサット画像８０に回転を適用する前および適用の後の、我々の方法および以前の方法によって推定された変位ベクトルの例を示す。図１２ｂは、数学的に計算され、提示される理想的な回転を示す。我々の方法によって推定された回転パターン１２２は、理想的な回転１２１に最も近いパターンである。この画像の平均レジストレーション誤差は、我々の方法に関して０．０８３４画素であった。テーブル４は、０．５°から３°までの回転値に関する、既存の方法および説明される方法のＲＥＲを要約するものである。我々の提案した方法は、すべての回転に関して最小の平均ＲＥＲを有する。他の方法は、ある度数を越える回転の変位を見つけることができない。これは、部分的には、１°を超える回転が、２画素を超える変位を引き起こすからである。方法は、大きい変位を頑健には推定しない場合に、この回転を正確に計算しまたは解決することができなくなる。
【０２００】
テーブル４
回転に関する２４００個すべてのランドサット部分画像（１２８画素×１２８画素）の平均レジストレーション誤差（画素単位）（ＲＥＲ）。最後の行は、我々の提案される方法の平均レジストレーション誤差に関して定量化された平均レジストレーション誤差の相対差を提供する。
【表３】
2019517073.tif 000018
【０２０１】
いくつかの場合に、特定の回転または変化が、ある方法に特に適する場合がある。ある方法が我々の提案した方法より小さい誤差を有した唯一の事例は、０．５°の回転の場合のＳｔｏｎｅの方法ＵＳ６６２８８４５であった。これの主な理由は、本方法の整数部分のＳＧＧＣ法とＣＣ法との間の差から生じる。たとえば、Ｓｔｏｎｅ他によって使用されるＣＣは、整数シフト（すなわち、１画素以上のシフト値）を見つけることに関して弱い。平均１．２８画素のシフトの場合に、シフトは、部分画像の多くで１画素未満であった。したがって、頑健な整数シフトの欠如は、小さい回転など、ある種の事例に関してＳｔｏｎｅの方法（ＣＣを使用し、画像のいくつかに関してシフトを記録しなかった場合がある）を妨げない。しかし、この大きさ範囲の外では、精度が失われる。特定の実施形態では、説明される方法は、この目的のために特に設計されていない場合であっても、画像をより小さい部分画像に分割することによって回転を推定することができる。
【０２０２】
ガウス雑音の存在下で、整数誤差メトリックおよび副画素誤差メトリックが、画像回転試験のために評価され得る。このために、異なるレベルの画像特徴を有する２つのランドサット画像（すなわち、パリ（図５（ａ））およびミシシッピ（図５（ｂ））が選択された。パリのランドサット画像は、ミシシッピのランドサット画像と比較して、より多くの特徴を有する。特徴の密度（および対応する特徴サイズ）は、本方法の精度に影響する可能性がある。０．５°と３°との間の回転が、各画像に適用され、我々の方法の変位レジストレーション誤差および２つの誤差メトリックが、各回転角で計算された。その結果は、それぞれパリおよびミシシッピのランドサット画像に関してテーブル５およびテーブル６に要約されている。その結果は、本方法の整数部分が、ガウス雑音よりも回転に敏感であることを示す。その結果は、両方の画像に関しておおまかに同様の誤差メトリックを示し、異なるテクスチャを有する異なる画像に対する本システムの頑健性を証明する。
【０２０３】
テーブル５
パリのランドサット画像に関する我々のアルゴリズムの整数部分および副画素部分の誤差メトリック
【表４】
2019517073.tif 000019
【０２０４】
テーブル６
ミシシッピのランドサット画像に関する我々のアルゴリズムの整数部分および副画素部分の誤差メトリック
【表５】
2019517073.tif 000020
【０２０５】
本発明の実施形態は、反復最適化プロセスを使用する方法またはアップサンプリング法と比較してかなり低い計算コストを有することができる。本発明の実施形態は、特異値分解および最適化法の使用より少ない計算を必要とする。ＳＶＤおよび２Ｄ  ＦＦＴの計算複雑さは、それぞれＯ（Ｎ＾３）および（Ｎ＾２ｌｏｇＮ）である。これは、１２４画素に関して６０．７４倍高速の計算を示唆する。我々の提案される方法の相対的に低い計算コストおよびその並列の性質は、この方法を、近リアルタイム応用例のためのＧＰＵおよびＦＰＧＡでのハードウェア実施に適するものにする。
【０２０６】
本発明をハードウェア実施に適するものにする本発明の態様は、
低い計算複雑さ（本方法は複雑な演算を全く含まない）と、
最適化プロセスなし（反復最適化は、通常、ハードウェア・アクセラレータでの実施が難しい）と、
低いメモリ・ストレージ要件（我々の方法は、小さいメモリ空間を要求する）と
を含む。並列性は、ハードウェア（たとえば、マルチコア・プロセッサ、ＦＰＧＡ、ＧＰＵ）で処理速度を向上させる（すなわち、計算時間を減少させる）最良の形の１つである。並列の性質を有する方法は、完全にまたは少なくとも実質的に独立な、分離可能な部分（並列に走行することができる）からなる。説明される方法は、レジストレーションが各部分画像内で独立なので、ある画像の部分画像に同時に並列に適用され得る。この並列性は、並列ハードウェア・アクセラレータ（たとえば、ＦＰＧＡ、ＧＰＵ、およびマルチコア・プロセッサ）でのかなりの速度向上につながる。
【０２０７】
並進シフト試験に加えて、回転の存在下で変位を測定する副画素画像レジストレーション方法を評価する方法が説明された。本方法は、回転パターンおよび変位を見つけることにおいて既存の方法より性能が優れている。
【０２０８】
本発明の実施形態は、副画素並進シフトの測定用に設計され、その高いレベルの精度および頑健性のゆえに、非剛体変換に関する粗から微細までの画像レジストレーション技法に組み込まれる能力を有する。
【０２０９】
本発明の特定の実施形態では、画像レジストレーションは、様々な柔らかい材料および固い材料の表面変形測定に適用され得る。たとえば、軟組織（たとえば、皮膚および胸）の表面変形は、直接測定によってまたはカメラ画像を使用して計算することが難しい。しかし、上で説明した画像レジストレーション方法の実施形態の使用は、リモート画像からの正確な測定を可能にする。一実施形態では、柔らかい材料または固い材料（たとえば、軟組織）の一連の写真が撮影され得、画像レジストレーション方法が、変位を見つけるために画像のそれぞれを正確に整列させるのに使用される。その後、柔らかい材料または固い材料（たとえば、軟組織）の変化または変形が、画像内の目立つ変化から検出され得る。本明細書では軟組織およびカメラ画像に関連して説明されるが、本技法は、様々な画像モーダリティ（たとえば、ＭＲＩ、ＣＴスキャン、超音波、顕微鏡）を使用して他の器官（たとえば、心臓組織、脳、および肝臓）内の変化、変形、または動きを見つける他の医療応用例に同様に適用され得る。いくつかの実施態様では、画像は、可視スペクトルまたは光学スペクトルに含まれない場合がある。たとえば、ＭＲＩ画像およびＸ線画像が、レジスタされ得る。本方法は、たとえば患者の運動または呼吸によって引き起こされる動き補償を可能にするために医療画像に適用され得る。
【０２１０】
さらなる実施形態では、本画像レジストレーション技法が、産業イメージング・システムに適用される。たとえば、本画像レジストレーション技法は、果物などの食品を分類し、等級付けする機械類に適用され得る。この技法は、システムを高フレーム・レート・カメラとインターフェースすることを必要とする。別の産業の例は、弾性材料のポアソン比の測定である。本システムの実施形態は、低い表面特徴が存在する場合に、変化に基づく頑健な処理のために適合され得る。本発明の実施形態は、高いまたは低い表面特徴の画像に関する要件に適用するために変更され得る。さらなる例では、本画像レジストレーション方法は、流量測定、具体的には粒子画像速度測定（ＰＩＶ）ビデオ（画像）に適用され得る。流量測定は、画像内の変化を見ることによって、デバイスのガイドまたはチャネルを通る流体の運動を追跡することができる。様々な応用例が望まれるが、これらは、本発明の範囲を限定しない。さらに、本画像レジストレーション方法は、精度、速度、および複雑さをトレード・オフすることによって特定の実施形態に合わせて調整され得る。
【０２１１】
図１４に示された方法は、複素部分画像の使用が、要求されるフーリエ変換の個数の低減を可能にするので（フーリエ変換は相対的に計算的に高価なので）、特に２次元画像と共に良好に働く。本方法は、レジストレーション目的関数または一致判断基準が、本方法が使用されるたびに同一の形で形成されなければならないことも仮定する。すなわち、レジストレーション一致判断基準は、部分画像の一次導関数の和の相互相関または強度交差微分などの様々な他の画像特徴として形成されなければならない。適切な平滑化関数の選択は、上では、雑音の低減と真の基礎になる値を保存するために適当な周波数プロパティを有することとのバランスとして説明された。適切な平滑化関数は、その平滑化関数のサポートの関数とすることもできる。サポートは、通常は、たとえば強度勾配（または別の画像特徴）を見つけるために選択された点の個数に直接に関係する。１点の近傍が、小さいサポートの例である。相対的に小さいサポートは、畳み込み演算の計算コスト（サポートのサイズに正比例する）を低減するが、適当な周波数応答の選択に関する制限要因である。しかし、説明される方法では、我々は、周波数領域でフィルタを適用することによって、大きいサポートを有するフィルタの計算コストの問題に対処した。他の実施形態では、より複雑な平滑化関数が要求される場合があり、あるいは、画像または画像特徴からのより高次元の入力が使用される場合がある。
【０２１２】
図１７の実施形態では、Ｎ−Ｄフィルタ（１７１）が、周波数領域平滑化微分器カーネル（１７２）の形で提示される。
【０２１３】
図２１ａは、３Ｄ試験画像（体積）内の相互相関（ＣＣ）と比較したフィルタリングされた相互相関（ＦＣＣ）を用いる方法の実施形態の利点を示す。これらの測定に使用されたＮ−Ｄフィルタは、３Ｄ  Ｓａｖｉｔｚｋｙ−Ｇｏｌａｙフィルタ（１７１）であり、逆ＦＦＴ（１７４）の出力は、等しく１の重みを付けられた。図２１ｂは、３Ｄ試験体積内の相互相関（ＣＣ）と比較したフィルタリングされた相互相関（ＦＣＣ）の雑音頑健性を示す。図２１ｃは、様々なサイズを有する試験体積について説明された方法の実施形態の精度を示す  は下のプロットに示される。
【０２１４】
図２２は、画像レジストレーション方法の概略図を示す。画像レジストレーション方法１９０は、好ましくは、プロセッサ１９１上または関連するメモリ１９２もしくはストレージ１９３内に含まれる。プロセッサは、画像ソースからを含む複数の入力を受け取ることができる。画像ソースは、２Ｄおよび／または３Ｄのカメラまたは測定デバイスを含む。これらは、Ｘ線１９５、３Ｄスキャナ１９６、またはカメラ１９７などの患者測定技法を含む。しかし、本方法は、これらの技法または画像ソースに限定されない。いくつかの実施形態では、画像ソースは、カメラ１９７に関して示されているように、異なる時刻の同一のソースとすることができる。本画像レジストレーション方法は、ユーザが、たとえば特徴抽出器フィルタを選択することによって、画像レジストレーション方法を適合させまたは制御することを可能にするユーザ・インターフェース１９８をも有することができる。本画像レジストレーション方法は、本画像レジストレーション方法を実行する前に画像を記憶するのにメモリ１９２またはストレージ１９３を使用することができる。レジスタされた画像、これらの組合せ、またはレジスタされた画像の組合せからの出力は、モニタまたはディスプレイ手段１９４上に表示され得る。
【０２１５】
当業者は、説明された方法またはその一部が、ＬＡＮ、ＷＡＮ、および／またはインターネットを含む有線および／または無線のネットワーク内で通信可能に結合された、デスクトップ・パーソナル・コンピュータもしくはラップトップ・パーソナル・コンピュータ、モバイル・デバイス、サーバ、および／またはその組合せなどの汎用デジタル・コンピューティング・デバイスによって実行されることが意図されていることを了解するであろう。具体的には、本システムは、好ましくは、汎用コンピュータ内のプロセッサを使用して実施されるが、代替アーキテクチャが、本発明の範囲から逸脱せずに可能である。説明されるプロセッサ手段またはプロセッサ１９０は、ＧＰＵ、ＦＰＧＡ、デジタル信号プロセッサ（ＤＳＰ）、マイクロプロセッサ、または他の処理手段とすることができる。本システムは、レジストレーションのために画像を記録するカメラまたは他の画像記録デバイス１９５、１９６、１９７をも含むことができる。そのようなデバイスまたは画像は、ＭＲＩ画像、超音波、ＣＴスキャン、顕微鏡画像、および衛星画像を含む。他の実施形態では、本画像レジストレーション方法は、カメラ内またはカメラに関連するプロセッサ内に記憶され得る。
【０２１６】
本発明の方法を実施するプログラム・ソフトウェアからの命令に従って特定の機能を実行するようにプログラムされた後に、そのようなデジタル・コンピュータ・システムは、事実上、特に本発明の方法に対する特殊目的コンピュータになる。このために必要な技法は、コンピュータ・システムの当業者に周知である。
【０２１７】
前述から、整数シフトおよび副画素シフトに関する改善された精度を提供する、画像レジストレーションの方法およびシステムが提供されることがわかる。
【０２１８】
文脈がそうではないことを明瞭に要求しない限り、この説明全体を通じて、単語「含む（ｃｏｍｐｒｉｓｅ、ｃｏｍｐｒｉｓｉｎｇ）」および類似物は、排他的または網羅的な意味ではなく包含的な意味で、すなわち、「〜を含むがこれに限定されない」という意味で解釈されなければならない
【０２１９】
本発明が、例によって、その可能な実施形態を参照して説明されたが、本発明の範囲から逸脱せずに、変更または改善が本発明に対して行われ得ることを理解されたい。本発明は、部分、要素、または特徴のうちの２つ以上の任意のまたはすべての組合せで、個別にまたは集合的に、本願の明細書で参照されまたは示された前記部分、要素、または特徴に存すると広義に言われることも可能である。さらに、既知の同等物を有する本発明の特定の構成要素または整数に対して言及が行われる場合に、そのような同等物は、個別に示されたかのように本明細書に組み込まれる。
【０２２０】
本明細書全体を通じての既存の方法または従来技術のすべての議論は、決して、そのような従来技術が広く知られているか、当分野における一般的な全般的知識の一部を形成することの容認と解釈されてはならない。
【図１】
2019517073.tif 000021
【図２】
2019517073.tif 000022
【図３】
2019517073.tif 000023
【図４】
2019517073.tif 000024
【図５】
2019517073.tif 000025
【図６】
2019517073.tif 000026
【図７】
2019517073.tif 000027
【図８】
2019517073.tif 000028
【図９】
2019517073.tif 000029
【図１０】
2019517073.tif 000030
【図１１】
2019517073.tif 000031
【図１２】
2019517073.tif 000032
【図１３】
2019517073.tif 000033
【図１４】
2019517073.tif 000034
【図１５】
2019517073.tif 000035
【図１６】
2019517073.tif 000036
【図１７】
2019517073.tif 000037
【図１８】
2019517073.tif 000038
【図１９】
2019517073.tif 000039
【図２０】
2019517073.tif 000040
【図２１】
2019517073.tif 000041
【図２２】
2019517073.tif 000042
【国際調査報告】
2019517073.tif 000043
2019517073.tif 000044
2019517073.tif 000045