(57)【要約】【課題】画像の対を位置合わせする方法と、プログラムを具現化してこの方法を実施するプログラム記憶装置とを提供して、例えばマルチモダリティで得られた３次元医用画像を精確にレジストレーションすること。【解決手段】画像の対を位置合わせする方法において、第1画像および第２画像を有する画像の対を用意し、ここで該画像には、３次元空間におけるピクセルのドメインに相応する複数の輝度が含まれており、前記の第１画像および第２画像の両方にて顕著特徴領域を識別し、ここで各領域には空間スケールが関連付けられており、各領域の中心点によって特徴領域を表し、局所的な輝度に基づいて、一方の画像の特徴点と、他方の画像の特徴点とをレジストレーションし、前記特徴対を類似性の尺度で順序付け、前記の中心点をサブピクセル精度に改善することによって特徴対の結合形対応付け集合を最適化する。【選択図】図６
2007054636.tif 000002
【特許請求の範囲】
【請求項１】
画像の対を位置合わせする方法において、
第1画像および第２画像を有する画像の対を用意し、ここで該画像には、３次元空間におけるピクセルのドメインに相応する複数の輝度が含まれており、
前記の第１画像および第２画像の両方にて顕著特徴領域を識別し、ここで各領域には空間スケールが関連付けられており、
各領域の中心点によって特徴領域を表し、
局所的な輝度に基づいて、一方の画像の特徴点と、他方の画像の特徴点とをレジストレーションし、
前記特徴対を類似性の尺度で順序付け、
前記の中心点をサブピクセル精度に改善することによって特徴対の結合形対応付け集合を最適化することを特徴とする、
画像の対を位置合わせする方法。
【請求項２】
さらに、
一方の画像に対する顕著特徴領域中心点をｋＤツリーで表し、
特徴毎に当該ｋＤツリーを調べて最隣接特徴の集合を探し、
顕著性値が低く、各特徴のスケール内に中心点を有する最隣接特徴を当該ツリーから除去して、前記の画像にて顕著特徴領域の実質的な均一な分布を得る、
請求項１に記載の方法。
【請求項３】
前記の空間スケールは、特徴領域を含む球の半径である、
請求項１に記載の方法。
【請求項４】
前記のｋＤツリーは、顕著特徴領域中心点の画像ピクセルインデックスをリーフとして使用し、
ここで特徴領域と、最隣接特徴領域との距離は、画像インデックスユニットである、
請求項２に記載の方法。
【請求項５】
局所的な輝度に基づく特徴点のレジストレーションではさらに、
前記の第１画像と第２画像との間のＩＣＰ変換を使用して初期レジストレーションを推定し、
第２画像のすべての特徴を第１画像の座標空間に変換し、
変換した当該の特徴をｋＤツリーに格納し、
上記の第１画像の特徴毎にこのツリーを調べてあらかじめ定めた選択判定基準に基づき第２画像にて最隣接特徴を選択し、
第１画像および第２画像の特徴からなる当該の選択された特徴対の並進不変性、回転不変性および大域的な画像類似性尺度をテストして、
大域的な画像類似性尺度値にしたがって当該の選択された特徴対を順序付ける、
請求項１に記載の方法。
【請求項６】
前記のＩＣＰ変換により、特徴点からなる各集合間の平均誤差を最小化する、
請求項５に記載の方法。
【請求項７】
前記の並進不変性のテストには、
2007054636.tif 000037
を推定することが含まれ、
ここでｐiおよびｐjは、物理空間におけるｉ番目の第１画像特徴およびｊ番目の第２画像特徴の中心座標である、
請求項５に記載の方法。
【請求項８】
前記の回転不変性のテストには
【数１】
2007054636.tif 000038
を推定することが含まれ、
ここで（ｆi，ｆj）はそれぞれ前記の第１画像および第２画像における特徴対を表しており、また
【数２】
2007054636.tif 000039
であり、ここでＨは、
【数３】
2007054636.tif 000040
と定められる、空間スケールｓを有するボクセル位置ｘの周りの球状隣接特徴領域（spherical neighborhood feature region）ｆs内の、画像輝度値を基準にしたエントロピーを表し、
ただしｐ(ｉ，ｓ，ｘ)は、ｆに含まれる画像輝度値ｉに対する確率密度関数であり、
Ｈ(ｆi，ｆj)は、
【数４】
2007054636.tif 000041
と定められる同時微分エントロピーであり、
ここでｐ(ｆi，ｆj)は、特徴領域ｆiおよびｆjにおける画像輝度の同時確率密度であり、ＩおよびＪはそれぞれ第１画像および第２画像において発生し得る輝度値集合にて値をとる、
請求項５に記載の方法。
【請求項９】
前記の大域的な画像類似性尺度のテストには、
【数５】
2007054636.tif 000042
を推定することが含まれ、
ここでＩrは第１画像を表し、
2007054636.tif 000043
は、第１画像の座標空間への第２画像の変換を表し、また
【数６】
2007054636.tif 000044
であり、ここでＨは、
【数７】
2007054636.tif 000045
と定められる、空間スケールｓを有するボクセル位置ｘの周りの、前記の画像のうちの１画像内の画像輝度値を基準にしたエントロピーを表し、
ただしｐ(ｉ，ｓ，ｘ)は、Ｉに含まれる画像輝度値ｉに対する確率密度関数であり、
2007054636.tif 000046
は
【数８】
2007054636.tif 000047
と定められる同時微分エントロピーであり、ここで
2007054636.tif 000048
は、Ｉrおよび
2007054636.tif 000049
における画像輝度の同時確率密度であり、ＩおよびＪはそれぞれ第１画像および第２画像において発生し得る輝度値集合において値をとり、これによってＬglobalが、第１画像および第２画像の重なりドメイン全体にて評価される、
請求項５に記載の方法。
【請求項１０】
特徴対の結合形対応付け集合の最適化ではさらに、
前記の類似性尺度にしたがい、最も類似している特徴対によって前記の結合形対応付け集合を初期化し、
結合形対応付け集合と、該結合形対応付け集合にまだ含まれていない各特徴対との和集合に対して類似性尺度を推定し、
当該の和集合の類似性尺度を最大にする特徴対を選択して、
最大値を与える特徴対と、結合形対応付け集合との和集合の類似性尺度が、この結合形対応付け集合の類似性尺度よりも大きい場合、前記の最大値を与える特徴対を、局所剛体変換によってサブピクセル精度でレジストレーションして、前記の結合形対応付け集合に加える、
請求項１に記載の方法。
【請求項１１】
ＩＣＰプロシージャを使用して前記の類似性尺度を最大化して、特徴対間のレジストレーション変換を計算する、
請求項１０に記載の方法。
【請求項１２】
最大値を与える特徴対と、結合形対応付け集合との和集合の類似性尺度が、当該結合形対応付け集合の類似性尺度以下である場合、類似性尺度を最大にする特徴対間のレジストレーション変換から計算したレジストレーション変換を供給する、
請求項１０に記載の方法。
【請求項１３】
画像の対を位置合わせする方法において、
第1画像および第２画像を有する画像の対を用意し、ここで該画像には、ピクセルからなる３次元空間の領域に相応する複数の輝度が含まれており、
前記の第１画像および第２画像の両方にて顕著特徴領域を識別し、ここで各領域には空間スケールが関連付けられており、
第１画像との第２画像との間でＩＣＰプロシージャを使用して初期レジストレーションを推定し、
第２画像のすべての特徴を第１画像の座標空間に変換し、
変換された特徴をｋＤツリーに格納し、第１画像の特徴毎に該ツリーを照会し、あらかじめ定めた選択判定条件に基づいて第２画像における最隣接特徴を選択し、
第１画像および第２画像の特徴からなる当該の選択した特徴対の並進不変性、回転不変性および大域的な画像類似性尺度をテストし、
当該の選択した特徴を、その大域的な画像類似性尺度値によって順序付けることを特徴とする、
画像の対を位置合わせする方法。
【請求項１４】
さらに、
各領域の中心点によって特徴領域を表し、
一方の画像に対する特徴領域中心点をｋＤツリーに格納し、
特徴毎に当該ｋＤツリーを調べて最隣接特徴の集合を探し、
顕著性値がより低く、各特徴のスケール内に中心点を有する最隣接特徴を当該ツリーから除去して、前記画像にて顕著特徴領域の実質的に均一な分布を得る、
請求項１３に記載の方法。
【請求項１５】
さらに、
前記の類似性尺度にしたがい、最も類似している特徴対によって結合形対応付け集合を初期化し、
当該の結合形対応付け集合と、当該結合形対応付け集合にまだ含まれていない各特徴対との和集合に対して類似性尺度を推定し、
当該の和集合の類似性尺度を最大化する特徴対を選択して、
前記の最大値を与える特徴対と、結合形対応付け集合との和集合の類似性尺度が、当該結合形対応付け集合の類似性尺度よりも大きい場合、当該の最大値を与える特徴対を、局所剛体変換によってサブピクセル精度でレジストレーションして、当該特徴対を前記結合形対応付け集合に加え、
ここで大域的な画像類似性尺度は、
【数９】
2007054636.tif 000050
と定義され、ここでＩrは第１画像を表し、
2007054636.tif 000051
は前記の第１画像の座標空間への第２画像の変換を表し、また
【数１０】
2007054636.tif 000052
であり、ここでＨは、
【数１１】
2007054636.tif 000053
と定められる、空間スケールｓを有するボクセル位置ｘの周りの一方の画像内の画像輝度値を基準にしたエントロピーを表し、
ただしｐ(ｉ，ｓ，ｘ)は、Ｉに含まれる画像輝度値ｉに対する確率密度関数であり、
2007054636.tif 000054
は、
【数１２】
2007054636.tif 000055
と定められる同時微分エントロピーであり、ここで
2007054636.tif 000056
画像Ｉrおよび
2007054636.tif 000057
における画像輝度の同時確率密度であり、ＩおよびＪはそれぞれ、第１および第２画像における発生し得る輝度値集合にて値をとり、これによってＬglobalは、第１および第２画像の重なりドメイン全体で評価される、
請求項１３に記載の方法。
【請求項１６】
コンピュータによって実行可能な命令からなるプログラムを具現化し、画像の対を位置合わせする方法ステップを実施する、コンピュータによって読み出し可能なプログラム記憶装置において、
前記の方法では、
第1画像および第２画像を有する画像の対を用意し、ここで該画像には、３次元空間におけるピクセルのドメインに相応する複数の輝度が含まれており、
前記の第１画像および第２画像の両方にて顕著特徴領域を識別し、ここで各領域には空間スケールが関連付けられており、
各領域の中心点によって特徴領域を表し、
局所的な輝度に基づいて、一方の画像の特徴点と、他方の画像の特徴点とをレジストレーションし、
前記の特徴対を類似性の尺度で順序付け、
前記の中心点をサブピクセル精度に改善することによって特徴対の結合形対応付け集合を最適化することを特徴とする、
コンピュータによって読み出し可能なプログラム記憶装置。
【請求項１７】
前記方法ではさらに、
一方の画像に対する顕著特徴領域中心点をｋＤツリーで表し、
特徴毎に当該ｋＤツリーを調べて最隣接特徴の集合を探し、
顕著性値がより低く、各特徴のスケール内に中心点を有する最隣接特徴を当該ツリーから除去して、前記の画像にて顕著特徴領域の実質的な均一な分布を得る、
請求項１６に記載のコンピュータによって読み出し可能なプログラム記憶装置。
【請求項１８】
前記の空間スケールは、特徴領域を含む球の半径である、
請求項１６に記載のコンピュータによって読み出し可能なプログラム記憶装置。
【請求項１９】
前記のｋＤツリーは、顕著特徴領域中心点の画像ピクセルインデックスをリーフとして使用し、
ここで特徴領域と、最隣接特徴領域との距離は、画像インデックスユニットである、
請求項１７に記載の、コンピュータによって読み出し可能なプログラム記憶装置。
【請求項２０】
局所的な輝度に基づく特徴点のレジストレーションではさらに、
前記の第１画像と第２画像との間のＩＣＰ変換を使用して初期レジストレーションを推定し、
第２画像のすべての特徴を第２画像の座標空間に変換し、
変換した当該の特徴をｋＤツリーに格納し、
上記の第１画像の特徴毎にこのツリーを調べてあらかじめ定めた選択判定基準に基づき第２画像にてて最隣接特徴を選択し、
第１画像および第２画像の特徴からなる選択された特徴対の並進不変性、回転不変性および大域的な画像類似性尺度をテストして、
大域的な画像類似性尺度にしたがって当該の選択された特徴対を順序付ける、
請求項１６に記載の、コンピュータによって読み出し可能なプログラム記憶装置。
【請求項２１】
前記のＩＣＰ変換により、特徴点からなる各集合間の平均誤差を最小化する、
請求項２０に記載の、コンピュータによって読み出し可能なプログラム記憶装置。
【請求項２２】
前記の並進不変性のテストには、
2007054636.tif 000058
を推定することが含まれ、
ここでｐiおよびｐjは、物理空間におけるｉ番目の第１画像特徴およびｊ番目の第２画像特徴の中心座標である、
請求項２０に記載の、コンピュータによって読み出し可能なプログラム記憶装置。
【請求項２３】
前記の回転不変性のテストには
【数１３】
2007054636.tif 000059
を推定することが含まれ、
ここで（ｆi，ｆj）はそれぞれ前記の第１画像および第２画像における特徴対を表しており、また
【数１４】
2007054636.tif 000060
であり、ここでＨは、
【数１５】
2007054636.tif 000061
と定められる、空間スケールｓを有するボクセル位置ｘの周りの球状隣接特徴領域ｆs内の、画像輝度値を基準にしたエントロピーを表し、
ただしｐ(ｉ，ｓ，ｘ)は、ｆに含まれる画像輝度値ｉに対する確率密度関数であり、
Ｈ(ｆi，ｆj)は、
【数１６】
2007054636.tif 000062
と定められる同時微分エントロピーであり、
ここでｐ(ｆi，ｆj)は、特徴領域ｆiおよびｆjにおける画像輝度の同時確率密度であり、ＩおよびＪはそれぞれ第１画像および第２画像において発生し得る輝度値集合にて値をとる、
請求項２０に記載の、コンピュータによって読み出し可能なプログラム記憶装置。
【請求項２４】
前記の大域的な画像類似性尺度のテストには、
【数１７】
2007054636.tif 000063
を推定することが含まれ、
ここでＩrは第１画像を表し、
2007054636.tif 000064
は、第１画像の座標空間への第２画像の変換を表し、また
【数１８】
2007054636.tif 000065
であり、ここでＨは、
【数１９】
2007054636.tif 000066
と定められる、空間スケールｓを有するボクセル位置ｘの周りの、前記の画像のうちの１画像内の画像輝度値を基準にしたエントロピーを表し、
ただしｐ(ｉ，ｓ，ｘ)は、Ｉに含まれる画像輝度値ｉに対する確率密度関数であり、
2007054636.tif 000067
は
【数２０】
2007054636.tif 000068
と定められる同時微分エントロピーであり、ここで
2007054636.tif 000069
は、Ｉrおよび
2007054636.tif 000070
における画像輝度の同時確率密度であり、ＩおよびＪはそれぞれ第１画像および第２画像において発生し得る輝度値集合において値をとり、これによってＬglobalが、第１画像および第２画像の重なり領域全体にて評価される、
請求項２０に記載の、コンピュータによって読み出し可能なプログラム記憶装置。
【請求項２５】
特徴対の結合形対応付け集合の最適化ではさらに、
前記の類似性尺度にしたがい、最も類似している特徴対によって前記の結合形対応付け集合を初期化し、
結合形対応付け集合と、該結合形対応付け集合にまだ含まれていない各特徴対との和集合に対して類似性尺度を推定し、
当該の和集合の類似性尺度を最大にする特徴対を選択して、
最大値を与える特徴対と、結合形対応付け集合との和集合の類似性尺度が、この結合形対応付け集合の類似性尺度よりも大きい場合、前記の最大値を与える特徴対を、局所剛体変換によってサブピクセル精度でレジストレーションして、前記の結合形対応付け集合に加える、
請求項１６に記載の、コンピュータによって読み出し可能なプログラム記憶装置。
【請求項２６】
ＩＣＰプロシージャを使用して前記の類似性尺度を最大化して、特徴対間のレジストレーション変換を計算する、
請求項２５に記載の、コンピュータによって読み出し可能なプログラム記憶装置。
【請求項２７】
最大値を与える特徴対と、結合形対応付け集合との和集合の類似性尺度が、当該結合形対応付け集合の類似性尺度以下である場合、類似性尺度を最大にする特徴対間のレジストレーション変換から計算したレジストレーション変換を供給する、
請求項２５に記載の、コンピュータによって読み出し可能なプログラム記憶装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、３Ｄデジタル医用画像用のレジストレーション法に関するものであり、殊に画像の対を位置合わせする方法およびコンピュータによって実行可能な命令からなるプログラムを具現化してこの方法を実施するプログラム記憶装置に関する。また本発明は、２００５年８月２４日に提出されたXu等による"Method and System for Robust Salient Region Based Registration of 3D Medical Images"という名称の米国暫定特許第６０／７１０，８３４号に優先権を主張するものであり、その内容をここに文献として援用する。
【背景技術】
【０００２】
医用画像処理においてレジストレーションは、２つ以上の画像の空間的な位置の間のマッピングを得る基本的な作業になっており、またさまざまな応用に使用可能である。位置合わせ変換（alignment transformation）に対する主要な要求は、対応する画像内容を最適に重ね合わせることである。既存のレジストレーションアプローチは、特徴ベースト（例えばランドマーク、エッジ、マーカ）、輝度ベースト、またこれらの２つの手法の様相が組み合わせられたハイブリッドな手法に分類することができる。ランドマークベーストの手法は、特徴から離れた領域においてうまく位置合わせができないことがある。画像の輝度だけに基づくアプローチは、最適化処理の性質に起因して、理想解からかけ離れた局所的な最適解に陥ってしまうことがある。ハイブリッド法は、これらの複数の特性の組み合わせを使用しており、また画像輝度または幾何学的な情報だけでは測定に対する精確なベースが得られない、異なるモダリティからの画像をレジストレーションするのには殊に有利となり得る。例えば、領域ラベル付け情報（region labeling information）からなる付加的な情報チャネルと対になった相互情報量により、ＭＲおよびＰＥＴ画像のレジストレーション結果を改善することができる。ハイブリッドレジストレーション法は、例えば、コントラストの低い腹部、血漿ゲル電気泳動またはプロテインイメージングなどにおいて既知である。別の応用分野は、例えば、画像または対象体分析に有利な標準データベース（normative database）またはアトラスの作成、医師が疾患の発達についての知識を得られるようにする１患者におけるまたは患者間の検討、またはガン治療中の時間ベースのフォローアップによる検討などである。
【０００３】
同じ被験体に対して異なる画像システムを使用することによって、より一層多くの情報が得られるが、その一方で適切な解釈のためにマルチモダリティのレジストレーション法が必要になる。補助情報の追加は、さまざまな医用イメージングシステムによって行われ、これらはおおまかには２つの主要なカテゴリに分類される。すなわち、形態学的な情報を抽出する解剖学的イメージング（例えば、Ｘ線、コンピュータトモグラフィ（ＣＴ）、磁気共鳴イメージング（ＭＲＩ）、超音波（ＵＳ））と、下部の解剖学的構造の代謝についての情報を視覚化する機能的イメージング（例えば、単光子放射型コンピュータ断層撮影法（SPECT Single Photon Emission Computed Tomography）、陽電子放射断層撮影（ＰＥＴ Positron Emission Tomography）、機能的ＭＲＩ（fMRI））とがある。マルチモーダル画像レジストレーションにおいて、異なるタイプの画像を組み合わせることは医師にとって有利である。例えば、ＣＴ画像によって良好な空間解像度が得られるのに対して、ＰＥＴ画像によって下部組織の機能が示されるのである。したがってＣＴ画像における機能的な情報の欠如は、空間解像度が欠如している対応するＰＥＴ画像との融合によって補償することができるのである。
【発明の開示】
【発明が解決しようとする課題】
【０００４】
本発明の課題は、画像の対を位置合わせする方法と、コンピュータによって実行可能な命令からなるプログラムを具現化してこの方法を実施するプログラム記憶装置とを提供して、例えばマルチモダリティで得られた３次元医用画像を精確にレジストレーションすることである。
【課題を解決するための手段】
【０００５】
上記の画像対を位置合わせする方法についての課題は、本発明の請求項１により、画像の対を位置合わせする方法において、第1画像および第２画像を有する画像の対を用意し、ここでこれらの画像には、３次元空間におけるピクセルのドメインに相応する複数の輝度が含まれており、上記の第１画像および第２画像の両方にて顕著特徴領域を識別し、ここで各領域には空間スケールが関連付けられており、各領域の中心点によって特徴領域を表し、局所的な輝度に基づいて、一方の画像の特徴点と、他方の画像の特徴点とをレジストレーションし、上記の特徴対を類似性の尺度で順序付け、上記の中心点をサブピクセル精度に改善することによって特徴対の結合形対応付け集合を最適化することによって解決される。
【０００６】
また上記の画像対を位置合わせする方法についての課題は、本発明の請求項１３により、画像の対を位置合わせする方法において、第１画像および第２画像を有する画像の対を用意し、ここでこれらの画像には、３次元空間におけるピクセルのドメインに相応する複数の輝度が含まれており、上記の第１画像および第２画像の両方において顕著特徴領域を識別し、ここで各領域には空間スケールが関連付けられており、第１画像と第２画像との間でＩＣＰ（Iterative Closest Point）プロシージャを使用して初期レジストレーションを推定し、第２画像のすべての特徴を第１画像の座標空間に変換し、変換された特徴をｋＤツリーに格納し、第１画像の特徴毎にこのツリーを調べて、あらかじめ定めた選択判定条件に基づいて第２画像における最隣接特徴（nearest neighbor feature）を選択し、第１画像および第２画像の特徴からなるこれらの選択した特徴対の並進不変性、回転不変性および大域的な画像類似性尺度をテストし、これらの選択した特徴を、その大域的な画像類似性尺度値によって順序付けることを特徴とする、画像の対を位置合わせする方法によって解決される。
【０００７】
上記のコンピュータによって実行可能な命令からなるプログラムを具現化して上記の方法を実施するプログラム記憶装置についての課題は、本発明の請求項１６により、コンピュータによって実行可能な命令からなるプログラムを具現化し、画像の対を位置合わせする方法ステップを実施する、コンピュータによって読み出し可能なプログラム記憶装置において、上記の方法では、第1画像および第２画像を有する画像の対を用意し、ここでこれらの画像には、３次元空間におけるピクセルのドメインに相応する複数の輝度が含まれており、上記の第１画像および第２画像の両方において顕著特徴領域を識別し、ここで各領域には空間スケールが関連付けられており、各領域の中心点によって特徴領域を表し、局所的な輝度に基づいて、一方の画像の特徴点と、他方の画像の特徴点とをレジストレーションし、上記の特徴対を類似性の尺度で順序付け、上記の中心点をサブピクセル精度に改善することによって特徴対の結合形対応付け集合を最適化することを特徴とする、コンピュータによって読み出し可能なプログラム記憶装置によって解決される。
【発明を実施するための最良の形態】
【０００８】
ここで説明する本発明の実施例に一般的に含まれるのは、レジストレーションアプローチに対する方法およびシステムであり、このレジストレーションアプローチにより、２つの画像から３Ｄ領域が自動的に抽出され、これらの２つ画像の間で対応する対が探し出され、また剛体レジストレーション変換が行われて医用アプリケーションにおいて融合された結果が視覚化される。これらの本来の特徴のスケール、並進および回転について不変的な特性は３Ｄ用に適合されて、これによって下位の単一モダリティまたはマルチモダリティ３Ｄ医用画像間の変換が推定される。本発明の実施形態による方法により、機能ベースおよび輝度ベースの両方のアプローチの有利な側面が組み合わされ、また各画像において３Ｄ顕著領域特徴（３Ｄsalient region feature）の１集合が自動的に抽出され、対応付け集合（corresponence set）が推定され、アウトライアの除去を含めたサブピクセル精度の改善が行われる。ここでは領域成長ベースのアプローチ（region-growing based approach）が使用されて３Ｄ顕著領域特徴が抽出され、これにより、特徴のクラスタリングおよび相応する探索空間の複雑さが低減されるのである。
【０００９】
本発明の１実施形態によるアルゴリズムにより、ｋＤツリー構造を活用することによって顕著領域が高速にクラスタリングされ、また３Ｄ顕著特徴領域（３Ｄsalient feature region）の局所輝度ドリブンなレジストレーションが使用されて最適化が改善される。本発明の１実施形態による付加的な機能に含まれるのは、顕著領域の２つの対の重心にＩＣＰ（iterative closest point）アルゴリズムを使用して初期ポーズの推定を行うこと、および局所的な輝度ベースのレジストレーションを使用して重心を局所的に改善してサブピクセル精度を得ることである。
【００１０】
本発明の１実施形態によるアルゴリズムは、全自動、ロバストかつマルチモダリティ剛体画像レジストレーションアルゴリズムであり、これは３Ｄ画像を任意のポーズでレジストレーションすることができる。画像は、スケール不変のさまざまな３Ｄ特徴によってモデリングされ、幾何学的な構成上の制約を受ける。ここでは同じ画像における特徴の間での相対的な構成上の制約がフルに活用されて、複数の３Ｄ顕著特徴対の間の結合形対応付け（joint correspondence）が求められる。個々の特徴対の間の対応付けに比べ、結合形対応付けによって課せられる厳格な幾何学上の制約により、誤ったマッチングが格段に発生しにくくなる。結合形対応付け集合から推定される変換は、大域的な最適値に収束する。それは、新たな対を加えても大域的な画像位置合わせ品質をさらに改善することできないようになるまで、結合形対応付けがインクリメント式に加算されるからである。これは適切な収束判定条件によって達成される。
【００１１】
本発明の１様相によれば、画像の対を位置合わせする方法が提供され、この方法では、第１画像および第２画像を有する画像の対を用意し、ここでこれらの画像には、３次元空間におけるピクセルのドメインに相応する複数の輝度が含まれている。また第１画像および第２画像の両方において顕著特徴領域を識別し、ここで各領域には空間スケールが関連付けられている。さらに各領域の中心点によって特徴領域を表し、局所的な輝度に基づいて、一方の画像の特徴点と、他方の画像の特徴点とをレジストレーションし、この特徴対を類似性の尺度で順序付け、上記の中心点をサブピクセル精度に改善することによって上記の特徴対の結合形対応付け集合を最適化する。
【００１２】
本発明の別の実施形態によれば、上記の方法において一方の画像に対する顕著特徴領域中心点をｋＤツリーで表し、特徴毎にこのｋＤツリーを調べて最隣接特徴（nearest neighbor feature）の集合を探し、上記のツリーから、顕著性値がより低くまた各特徴のスケール内に中心点を有する最隣接特徴を取り除き、これによって上記の画像において顕著特徴領域の実質的な均一な分布を得る。
【００１３】
本発明のさらに別の１実施形態によれば、上記の空間スケールは、上記の特徴領域を含む球の半径である。
【００１４】
本発明のさらに別の１実施形態によれば、上記のｋＤツリーは、上記の顕著特徴領域中心点の画像ピクセルインデックスをリーフとして使用し、ここで特徴領域と、最隣接特徴領域との距離は、画像インデックスユニットで表される。
【００１５】
本発明のさらに別の１実施形態によれば、局所的な輝度に基づく特徴点のレジストレーションではさらに、上記の第１画像と第２画像との間のＩＣＰ変換を使用して初期レジストレーションを推定し、第２画像のすべての特徴を第１画像の座標空間に変換し、変換したこれらの特徴をｋＤツリーに格納し、上記の第１画像において特徴毎にこのツリーを調べ、あらかじめ定めた選択判定基準に基づいて第２画像において最隣接特徴を選択し、第１画像および第２画像の特徴からなるこの選択された特徴対の並進不変性、回転不変性、および大域的な画像類似性尺度をテストして、この選択された特徴対を大域的な画像類似性尺度値にしたがって順序付ける。
【００１６】
本発明のさらに別の１実施形態によれば、上記のＩＣＰ変換により、特徴点からなる各集合間の平均誤差が最小化される。
【００１７】
本発明のさらに別の１実施形態によれば、上記の並進不変性のテストには、
2007054636.tif 000003
を推定することが含まれ、ここでｐiおよびｐjは、物理空間におけるi番目の第１画像特徴およびj番目の第２画像特徴の中心座標である。
【００１８】
本発明の別の１実施形態によれば、上記の回転不変性のテストには
【数１】
2007054636.tif 000004
を推定することが含まれ、ここで（ｆi，ｆj）はそれぞれ上記の第１画像および第２画像における特徴対を表しており、また
【数２】
2007054636.tif 000005
であり、ここでＨは、
【数３】
2007054636.tif 000006
と定められる、空間スケールｓを有するボクセル位置ｘの周りの球状隣接特徴領域（spherical neighborhood feature region）ｆs内の、画像輝度値を基準にしたエントロピーを表し、ただしｐ(ｉ，ｓ，ｘ)は、ｆに含まれる画像輝度値ｉに対する確率密度関数であり、Ｈ(ｆi，ｆj)は、
【数４】
2007054636.tif 000007
と定められる同時微分エントロピー（joint differential entropy）であり、ここでｐ(ｆi，ｆj)は、特徴領域ｆiおよびｆjにおける画像輝度の同時確率密度（joint probability density）であり、ＩおよびＪはそれぞれ第１画像および第２画像において発生し得る輝度値集合において値をとる。
【００１９】
本発明の別の１実施形態によれば、上記の大域的な画像類似性尺度のテストには、
【数５】
2007054636.tif 000008
を推定することが含まれ、ここでＩrは第１画像を表し、
2007054636.tif 000009
は、第１画像の座標空間への第２画像の変換を表し、また
【数６】
2007054636.tif 000010
であり、ここでＨは、
【数７】
2007054636.tif 000011
と定められる、空間スケールｓを有するボクセル位置ｘの周りの、前記の画像のうちの１画像内の画像輝度値を基準にしたエントロピーを表し、ただしｐ(ｉ，ｓ，ｘ)は、Ｉに含まれる画像輝度値ｉに対する確率密度関数であり、
2007054636.tif 000012
は
【数８】
2007054636.tif 000013
と定められる同時微分エントロピーであり、ここで
2007054636.tif 000014
は、Ｉrおよび
2007054636.tif 000015
における画像輝度の同時確率密度であり、ＩおよびＪはそれぞれ第１画像および第２画像において発生し得る輝度値集合において値をとり、これによってＬglobalが、第１画像および第２画像の重なりドメイン全体で評価される。
【００２０】
本発明のさらに別の１実施形態によれば、特徴対の結合形対応付け集合の最適化ではさらに、上記の類似性尺度にしたがい、最も類似している特徴対によって上記の結合形対応付け集合を初期化し、結合形対応付け集合と、この結合形対応付け集合にまだ含まれていない各特徴対との和集合に対して類似性尺度を推定し、この和集合の類似性尺度を最大にする特徴対を選択して、最大値を与える特徴対と、結合形対応付け集合との和集合の類似性尺度が、この結合形対応付け集合の類似性尺度よりも大きい場合、上記の最大値を与える特徴対を、局所剛体変換（local rigid tranformation）によってサブピクセル精度でレジストレーションして、上記の結合形対応付け集合に加える。
【００２１】
本発明のさらに別の１実施形態によれば、ＩＣＰプロシージャを使用することによって上記の類似性尺度を最大化して、特徴対間のレジストレーション変換を計算する。
【００２２】
本発明のさらに別の１実施形態によれば、最大値を与える特徴対と、結合形対応付け集合との和集合の類似性尺度が、この結合形対応付け集合の類似性尺度以下である場合、類似性尺度を最大にする特徴対間のレジストレーション変換から計算したレジストレーション変換を供給する。
【００２３】
本発明の別の様相によれば、コンピュータに読み取り可能なプログラム記憶装置が提供され、ここでこのプログラム記憶装置は、コンピュータによって実行される命令のプログラムを具現化して上記の画像の対を位置合わせする方法の処理ステップを実施するものである。
【実施例】
【００２４】
ここで示す本発明の実施例に一般的に含まれるのは、スケール、回転および並進に対して不変である自動抽出した３Ｄ顕著領域特徴を構成的にマッチングする方法およびシステムである。したがって本発明にはさまざまな変形および択一的な形態が可能であるが、図面ではその特定の実施形態を例示して以下に詳しく説明する。しかしながら、本発明をここで示す特定の形態に制限する意図はなく、逆に本発明の精神および範囲内に入るすべての変形形態、同等の形態および択一的な形態が本発明に含まれることを理解されたい。
【００２５】
ここで使用される「画像」という語は、離散的な画像要素（例えば、２Ｄ画像に対するピクセル、３Ｄ画像に対するボクセル）からなる多次元データのことである。画像は、例えば、コンピュータトモグラフィ、磁気共鳴イメージング、超音波または当業者には公知の他の任意の医用イメージングシステムによって収集した被験体の医用画像とすることが可能である。画像は、医用でないコンテキスト、例えばリモートセンシングシステム、電子顕微鏡などから得ることも可能である。画像はＲ3からＲへの関数と考えることができるが、本発明の方法は、このような画像には限定されず、任意の次元の画像、例えば、２Ｄの像または３Ｄのボリュームに適用可能である。２次元または３次元の画像に対して、画像の領域はふつう２次元または３次元の矩形のアレイであり、各ピクセルまたはボクセルは、２つまたは３つの互いに垂直な軸の組を基準にしてアドレッシングすることができる。ここで使用する「デジタル」または「デジタイズされた」という語は、デジタル収集システムを介してまたはアナログ画像から変換を介して収集されたデジタルまたはデジタイズされた形式の画像またはボリュームのことを適宜表すこととする。
【００２６】
本発明の１実施形態によれば、顕著性の記述が使用されて、３Ｄ画像から特徴が自動的に抽出される。本発明の１実施形態による３Ｄレジストレーションアルゴリズムを要約するフローチャートが図６に示されている。これらのステップの詳細を以下に詳述する。図を参照すると、ステップ６０において球状領域を形成することによって３Ｄ顕著領域特徴が検出される。結果的に得られる各顕著領域特徴はつぎの特性を有する。すなわち、（１）領域の中心、（２）領域のスケール（半径）、および（３）領域の顕著性値である。
【００２７】
画像空間にわたって顕著領域特徴が均一に分布することは、ステップ６１において、一方の集合の特徴領域中心点をｋＤツリー構造に格納することによって達成される。これによって空間位置について高速に調べることができる。ステップ６２において、特徴毎にこのツリーを調べて最隣接特徴を探す。領域はそのスケール（半径）によって表されるため、顕著性がより低くかつ中心点がカレント特徴のスケール内にある特徴は、ステップ６３で除去される。結果的に得られる、クラスタリングされた顕著特徴の部分集合は、全体集合から取り除かれる。より大きな顕著性値を有し、顕著領域特徴のスケール内に含まれる顕著領域特徴はない。このアプローチによって、大きな顕著性値、すなわち局所的に最大の顕著性値を有し、画像領域に発生し得る顕著領域特徴をクラスタリングすることが回避される。したがってこれらの特徴に基づく後続の画像レジストレーションには偏りが発生して誤った解が得られる可能性がある。
【００２８】
回転不変性に対するパラメタは、３Ｄ顕著特徴領域の局所的な輝度ドリブンなレジストレーションによって推定される。局所変換のパラメタ探索空間は、回転パラメタに制限される。各組み合わせをテストすることによって相関していない２つの集合に含まれている顕著領域特徴間で対応付けを探すと、計算が極めて繁雑になる。２つの３Ｄ顕著領域特徴集合の場合、これらの特徴は、ステップ６４において中心点に還元（reduce）される。ステップ６５において結果的に得られる２つのポイントクラウド（point cloud）は、これらの集合間の平均二乗誤差を最小化するＩＣＰアルゴリズムによってレジストレーションされる。結果的得られた変換がステップ６６で使用されて、位置合わせ画像特徴集合（align image feature set）が基準画像の空間にマッピングされる。変換された位置合わせ特徴はステップ６７でｋＤツリー構造に格納され、これによって基準画像顕著領域特徴を基準にして、変換され最も近いＮ個の位置合わせ特徴を高速に調べることができる。初期のＩＣＰ変換によって、点集合間にロバストなマッピングが得られるとの前提に基づくと、仮定した対応付けに対するテストは、位置合わせ画像顕著特徴領域からなる全体集合の代わりに、基準特徴毎のＮ個の隣接に低減することができるのである。
【００２９】
顕著特徴領域の対応付け間でより一層高い精度を得るため、ステップ６８において特徴中心は、顕著特徴領域間の局所的な輝度ドリブンのレジストレーションによって位置合わせされる。この場合、期待値を最大化する最適化処理に対するパラメタ探索空間には、並進、回転およびスケールに対するパラメタが含まれる。
【００３０】
顕著領域特徴の基本な考え方は、所定のスケールについて信号の複雑性または局所的な予測不可能性が高いことを表すことあり、ここでスケールとはボクセルの周りの球状の領域の半径のことである。このアプローチにより、スケールの異なる円形領域に対し、シャノンのエントロピーに基づいて関心点が区別される。またこのアプローチでは、対応するが異なる解剖学的または機能的構造から得られるボクセルが、類似の顕著性値を有することが前提とされる。顕著性の記述の局所的な性質により、画像レジストレーションにとって大きな利点が得られる。すなわち、異なる画像間の対応する顕著領域特徴は、画像がオーバラップしない場合であっても全体的な空間変換に対して不変なのである。
【００３１】
SPECTおよびＰＥＴ画像の観察によって示されるのは、相応するＣＴ画像と比較すると、局所顕著性最大値の位置が、対応する関心構造内で局所的に並進し得ることである。これは局所剛体レジストレーション（local rigid registration）ステップによって扱うことができ、ここでこのステップには、画像内容を変形することなく、局所的な輝度相関に基づき、領域中心をサブピクセル精度で相応する位置に並進することが含まれる。このステップでは、前に述べたように、相応する特徴に対して顕著性値が類似しているといの基本的な前提が保持される。
【００３２】
顕著性は、画像輝度範囲Ｄに対してつぎのように定められる。すなわち、
ＡD(ｓp，ｘ) ＝ ＨD(ｓp，ｘ)・ＷD(ｓp，ｘ)
であり、ここでＨDは、スケール（半径）ｓを有するボクセル位置ｘの周りで球状隣接領域ＲS内の画像輝度値ｉ∈Ｄを基準にしたエントロピーを示し、
【数９】
2007054636.tif 000016
である。ここでｐ(ｉ，ｓ，ｘ)は、ＲSに含まれる画像輝度値に対する記述子ｉの確率密度関数（ＰＤＦ  probability density function）である。ＷD(ｓ，ｘ)は、スケールについてのＰＤＦ間の類似性の尺度である。これは、ＰＤＦの非類似性が増すのに伴って増大する。すなわち、
【数１０】
2007054636.tif 000017
である。
【００３３】
ｘにおけるＨDの局所的なピークになるスケールｓpは、
【数１１】
2007054636.tif 000018
によって得られる。
【００３４】
ボクセル毎に式（１）を解いた後、結果的に得られるのは、分析すべき入力画像と同じサイズの２つの一時画像である。すなわち、実際の顕著性値を含む画像と、式（２）から得られたスケール値を含む別の画像とが得られるのである。局所的に最適でありかつ最も記述的な顕著領域点は、領域成長サーチアプローチ（region-growing search approach）によって顕著性画像から抽出される。探索空間を低減するため、大域的顕著性閾値δを下限として使用する。本発明の１実施形態によれば、経験的に平均顕著性の半分に設定する、すなわち
2007054636.tif 000019
によって良好な結果が得られて重要でない領域を消去された。
【００３５】
領域特徴のｋＤツリー構造化に基づく最隣接点アルゴリズム（nearest point algorithm）により、局所的な顕著性最大値の位置を決定することができ、これによって顕著性値にしたがって順序付けられたボクセル位置のリストが得られる。このアプローチによって、局所的な最大値のクラスタリングが回避される。ここでこれは、例えば、大域的な閾値が適用されまた特徴が、低下する顕著性値だけにしたがって抽出される場合に発生し得るものである。ｋＤツリーは、抽出された局所顕著性最大値の領域中心のインデックスをリーフとして使用することによって生成される。特定の特徴についてのＫ個の最隣接特徴は、この特徴の領域中心インデックスについて調べることによって効率的に探すことができる。したがって返された特徴への距離は、物理的な距離ではなく、画像インデックス次元における距離である。スケールパラメタは、最小距離要求条件として使用可能である。すなわち問い合わせた特徴のスケール以下の距離を有しかつ返されたすへでの特徴および顕著性の低いすべての特徴は、特徴集合から取り除かれる。この制限は、集合全体に適用され、クラスタリングされた領域を除去することができる。結果的に得られる集合によって特定のサイズが要求される場合、距離の判定条件を満たす、顕著性のより低い特徴で上記のリストを埋めることができる。特徴は、その中心が、顕著性値のより高い特徴の領域内にない場合、集合に保持される。結果的に得られる３Ｄ顕著領域特徴の集合は、均一に分布し、また後続の特徴対応探索に対して良好な条件を有する初期集合が得られる。
【００３６】
図１（ａ）〜（ｅ）には、顕著特徴領域の集合の例が示されている。顕著特徴領域は図１（ａ）〜（ｂ）において白い円で示されている。図１（ａ）は、顕著領域特徴の集合にクラスタリングを行った効果を示しており、これに対して図１（ｂ）は、本発明の１実施形態によるアプローチによって選択された顕著領域特徴を示している。図１（ｃ），１（ｄ）および１（ｅ）はそれぞれＣＴ画像、ＰＥＴ画像およびＭＲ画像から抽出した後に視覚化した最も重要な３Ｄ領域特徴を示している。顕著領域特徴は、図１（ｃ）〜（ｅ）において球状のバブルとして示されている。これらのボリュームは、特定の伝達関数によってウィンドウィングされて、特徴の位置が３Ｄで視覚化される。これに対して抽出そのものは、すべての輝度範囲に対して実行される。
【００３７】
領域コンポーネントマッチングステップ（region component matching step）と称される３Ｄレジストレーションのつぎのステップにより、２つの画像の特徴間の、仮想的な対応の集合が推定される。Ｉrを基準画像、Ｉtをテンプレートとし、またＮrをＩrから抽出した特徴の数、ＮtをＩtから抽出した特徴の数とする。すべての仮想的な特徴の対応の集合は、Ｃ＝｛ｃｉ，ｊ｝であり、ただしｉ∈［１，Ｎr］，ｊ∈［１，Ｎt］，｜Ｃ｜＝Ｎr×Ｎtであり、またｃｉ，ｊ＝（ｆi，ｆj）は、Ｉrの特徴ｆiと、Ｉtの特徴ｆjとの対である。
【００３８】
パラメタ集合Θにより、２つの画像を位置合わせする変換Ｔが定められ、またこれは、ｆiおよびｆjの並進、スケールおよび回転について不変性に基づいて推定することができる。ｆiとｆjとの間の並進部分は、
2007054636.tif 000020
によって直接推定され、ただしｐiおよびｐjは、物理空間におけるｉ番目の基準特徴およびｊ番目のテンプレート特徴の中心位置である。スケールについての不変性はこの場合には不要である。それは本発明の実施形態による３Ｄ医用画像に対して、ボクセル寸法（voxel dimension）は、DICOM（Digital Imaging and Communications in Medecine）ヘッダ内で供給されるからである。回転不変性を得るため、輝度値に基づいて、３Ｄ顕著特徴領域の局所剛体レジストレーションによって回転パラメタが推定される。最適化は、回転パラメタ部分空間ΘRに制限され、また輝度の類似性の尺度、エントロピー相関係数によって行われ、これは正規化された相互情報量の特殊な形態、すなわち、
【数１２】
2007054636.tif 000021
であり、ここで同時微分エントロピーＨ（Ａ，Ｂ）は
【数１３】
2007054636.tif 000022
と定義できる。ここで積分のドメインはＲiおよびＲΘjにわたっており、ｐ（Ａ，Ｂ）は、領域ＡおよびＢにおける画像輝度の同時確率密度であり、ＩおよびＪはそれぞれ、ＩfおよびＩmにおける発生し得る輝度値集合において値をとる。この係数により、オーバラップドメインに対して改善したロバストネスが得られ、またいくつかの付加的な有利な特性が得られる。すなわち、画像間の依存性が増大することを示す値の増大およびその逆が得られるのである。したがって回転に対する不変性はΘRに対する最適化問題として定式化することができる。すなわち、
【数１４】
2007054636.tif 000023
である。
【００３９】
大域的な画像の類似性尺度Ｌglobalは、Ｍ個の各対の品質を推定するために使用される。すなわち
【数１５】
2007054636.tif 000024
であり、ここでＬglobalは、局所特徴領域だけではなく２つの画像の全体的なオーバラップドメインにおいて評価される。
【００４０】
大きな空間距離によって隔てられている特徴対は、対応するとは考えにくく、仮想的な対応付けの集合から取り除くことができる。これによって結合形対応付けに対する探索空間が低減される。したがって対応付け探索空間は、すべての対の組み合わせから、局所的に最も近い隣接特徴間だけの組み合わせに低減することができる。隣接集合は、基準およびテンプレート領域特徴集合間のＩＣＰ変換を計算することによって推定することができ、ここでこれは集合を、領域中心位置の周りのポイントクラウドと見なすことによって行われる。このＩＣＰアルゴリズムにより、特徴集合は、局所的な最小平均二乗誤差ＭＳＥ（minimal mean-square-error）に位置合わせされる。この結果は、すべてのテンプレート特徴を基準画像の座標空間に変換して、新たなｋＤツリーに格納するのに使用される。この場合、基準画像において顕著特徴毎に近似的に最も近い隣接特徴を、このツリーにおける高速探索によって決定することができる。各基準特徴に結合される、変換されたテンプレート特徴の最隣接特徴の数Ｎnは、この集合の全体濃度よりも格段に小さく、Ｎn≪Ｎtである。これにより、初期のＩＣＰ変換が実際の位置合わせ変換の良好な近似であるという前提、および特徴が比較的大きな距離で離れている場合、これらはあまり対応していないという前提に基づいて複雑さは、Ｎr×Ｎnに低減されるのである。Ｎr×Ｎnの特徴対の対応付けの品質がテストされ、ここでこれは並進不変量
2007054636.tif 000025
回転不変量
2007054636.tif 000026
大域的画像類似性尺度Ｌglobal(ｃｉ，ｊ)に基づいて行われ、これは位置合わせされた画像全体に局所的な特徴変換を適用することによって得られる。
【００４１】
本発明の実施形態にしたがって行われた実験において、Ｎn＝（１／１０）％Ｎtの隣接サイズが適用されて、結合形対応付けに対する初期探索空間を設けることができた。さらに仮想的な対応付けは、大域的類似度尺度Ｌglobalによって順次付けされる。これにより、推定された対応付け集合においてアウトライアはわずかになる。
【００４２】
図２（ａ）〜（ｄ）には、並進され、回転され、元々のスライスに重ね合わされたＣＴボリューム（左側部分）のスライスが示されている。円は、特定のスケールを有する顕著特徴領域を表している。順序付けがＬglobalにしたがって行われており、全体集合にはわずかなアウトライアしか含まれない。わかりやすくするため、アプローチ毎に最初の４つの対応付けだけが示されている。
【００４３】
このアルゴリズムの特徴対応付け探索において計算されるサイズＭの仮想的な対応付けの集合
2007054636.tif 000027
は、２つの画像間の変換Ｔを推定するために使用される。この変換の精度は十分ではない。それはそのパラメタは、離散的な画像グリッド位置に固定された特徴で計算されているからである。上述したように、精確に対応する空間位置に位置していない特徴対もある。したがって結果的に得られる集合は、変換をマイナスに偏らせるアウトライアおよび不確かさを含み得るのである。つぎのステップにおいてΘおよびＣは、サブピクセル精度の繰り返し処理において改善され、さらに精度の高い位置合わせが達成される。
【００４４】
サブピクセル精度で位置合わせされた特徴対を含みかつ理想的にはアウトライアを含まない結合形対応付け集合
2007054636.tif 000028
を最適化することが望ましい。最適化された結合形対応付け集合の要素は、ＩＣＰアルゴリズムに対する入力として使用されて、大域的な画像の類似性を最大化する変換が計算される。すなわち、
【数１６】
2007054636.tif 000029
である。
【００４５】
特徴対の数を少なく維持し、またレジストレーション効率を維持するため、ＥＭ（expectation-maximization）タイプアルゴリズムが、限られた数の繰り返しステップで使用される。変換
2007054636.tif 000030
は、各繰り返しにおいて徐々に改善される結合形対応付け集合Ｊk＿Ｊから計算される。Ｌglobalは、改善処理に対する収束条件として使用される。領域特徴対が一旦サブピクセル精度で局所的にレジストレーションされると、この特定の対の後続のレジストレーションによって、この対応付けの品質は向上されないため、これを無視することができる。したがって計算時間は、各繰り返しステップ中に繰り返して加えられる特徴対位置を改善するだけで節約することができる。
【００４６】
結合形対応付け集合およびレジストレーション変換の両方を最適化するＥＭタイプアルゴリズムのフローチャートは図３に示されている。このアルゴリズムは、Ｃの最上位の対、すなわち
2007054636.tif 000031
を含む結合形対応付け集合によって初期される。これらの初期の対応付けに対してふつう前のステップにおいて収集された特徴対の順序集合から最良の２つを使用することができる。局所剛体レジストレーションは、対応する顕著領域特徴をサブピクセル精度で改善するために使用される。図を参照すると、サブピクセルで改善された対応付けのカレント集合はステップ３１で得られる。ステップ３２の推定ステップでは、すべての
2007054636.tif 000032
に対して大域的類似性尺度Ｌglobal(Ｊ*∪ｃｉ，ｊ)を計算する。ステップ３３の最大化ステップでは、Ｌglobal(Ｊ*∪ｃｉ，ｊ)を最大化する
2007054636.tif 000033
を選択する。つぎにステップ３４では、最大のＬglobal(Ｊ*∪ｃｉ，ｊ)[Ｌglobal(Ｊ*)の場合、ステップ３５において変換ＴJ*が返され、このアルゴリズムは終了する。その他の場合、最大値を与える特徴対
2007054636.tif 000034
がステップ３５において局所剛体レジストレーションによってサブピクセル精度でレジストレーションされて特徴中心が改善される。改善された特徴対
2007054636.tif 000035
は、ステップ３７において集合Ｊ*に加えられ、すなわち、
2007054636.tif 000036
であり、変換ＴJ*がステップ３８で再計算される。ここでこれらのステップ３２〜３８は収束するまで繰り返される。
【００４７】
本発明の実施形態によるアルゴリズムをさまざまな一患者について３Ｄ医用画像についてテストした。異なる時間に収集した１１個のＰＥＴ−ＣＴボリューム対、治療の異なる段階における３つのＣＴ画像、およびハイブリッドスキャナから得られた１０個のSPECT−ＣＴボリューム対に対して測定を行った。本発明の１実施形態による方法は、さまざまなモダリティ、ノイズ、視野の変化、およびＰＥＴ−ＣＴ対のいくつかにおける画像輝度のアーチファクトなどを対処しなければならず、いくつかのスライスは、インポート中に補正されなかった異なる輝度スケールを有していた。図４（ａ）〜（ｃ）はそれぞれＰＥＴ−ＣＴ，ＣＴ−ＣＴおよびSPECT−ＣＴボリューム対に対して全体的な測定距離をｘ，ｙおよびｚ方向においてｃｍで、また標準偏差と共に示す表である。ＰＥＴ−ＣＴおよびＣＴ−ＣＴレジストレーション品質は、いくつかの関心対象点間の距離、すなわち右および左側の肺尖、心尖、肝臓の円形端部、左上および右上ならびに左下および右下の腎臓端部を医療の専門家が評価したものである。１０個のSPECT−ＣＴ画像は、最新式のハイブリッドスキャナによって収集したため、医師は、ｘ，ｙおよびｚ方向に１０〜５０ｍｍ変化させ、各軸の周りに５〜６０度回転することによってSPECT画像の手動で剛体的（rigidly）にレジストレーション解除（deregister）した。レジストレーションの後、いくつかの識別可能なランドマークが、ＣＴおよびSPECT画像において医療専門家によって選択された。
【００４８】
スライス間で輝度のスケーリングが変わることに起因するアーチファクトまたはノイズを含む実際の医用画像に対して複数の実験を行った。本発明の実施形態によるアルゴリズムをこれらのデータでテストするため、これらの問題はレジストレーションの前に処理しなかった。図５（ａ）〜（ｃ）は、本発明の１実施形態によるアルゴリズムを使用して収集したＰＥＴ−ＣＴ画像対、輝度アーチファクトを有するＣＴ−ＣＴ画像対、およびSPECT−ＣＴ画像対からそれぞれ得られ、融合式にレジストレーションされた結果画像からの３つのスライスを示している。後者のＣＴ画像は、限定された視野で収集され、大量のノイズを含むが、ここで提案したレジストレーションの精度は許容範囲になった。残りのレジストレーション誤りは、本発明の非剛性変換の実施形態によって処理することができる。
【００４９】
視覚化をおよび測定のための専用の視覚化ソフトウェアを使用して医療専門家によって結果を評価した。評価のため、医療専門家は、３Ｄの関心領域の重心を使用するのと、ランドマークから他のランドマークへの直接の位置と間で選択を行った。いくつかの付加的な測定ツールと共に融合式の視覚化を描画ソフトウェアに組み込むことによってこの作業を支援した。ＰＥＴ−ＣＴの場合、ｚ方向において標準偏差が高いことは明らかである。その理由は、収集モデル間の違いから発生し得る。ＣＴ画像は、１呼吸のスナップショットを示しているのに対して、ＰＥＴ画像は複数の呼吸周期にわたって収集されており、多かれ少なかれ平均的な呼吸動作を示している。横隔膜のこのような動作に起因して腹部のいつかの器官は持ち上げられたり押し下げられたりして、データサンプルに大きな偏差が見られる原因になっている。この実験に使用された本発明の１実施形態によるアルゴリズムは、剛体変換だけをモデリングし、このような局所的な変形をモデリングしない。ＣＴ−ＣＴデータに対しこのような作用はもはや支配的ではない。それは患者が２つの収集において理想的には同じように呼吸するからである。SPECT−ＣＴデータは、本質的に良好にマッチングし、SPECTにおけるユーザ定義の剛性変換によって局所的な変形は発生しない。したがってこれらのケースに対して良好なレジストレーションを期待することができる。
【００５０】
すべての結果において、所定の測定誤差が発生する。それは、医療専門家がさまざまなスライスビューにおいて複数の位置をクリックすることによって位置を手動で指定しなければならないからである。しかしながらこのタイプの評価についての実験を行うにあたって、いくつかの測定ステップにおいて関心対象点を指定する距離の平均差分は（１観察者においてまた観察者間において）３ｍｍを上回らない。
【００５１】
本発明は、ハードウェア、ソフトウェア、ファームウェア、特殊用途プロセッサまたはこれらの組み合わせなどのさまざまな形態で実現可能であることを理解されたい。１実施形態において本発明は、コンピュータ読み出し可能プログラム記憶装置上に具現化されるアプリケーションプログラムとしてソフトウェアでインプリメント可能である。このアプリケーションプログラムは、任意のアーキテクチャを含むマシンにアップロードして実行することができる。
【００５２】
図７は、本発明の１実施形態による３Ｄレジストレーション処理をインプリメントするコンピュータシステムの例のブロック図である。図７を参照すると、本発明を実施するコンピュータシステム７１は、殊に中央処理ユニット（ＣＰＵ）７２，記憶装置７３および入出力（Ｉ／Ｏ）インタフェース７４を含むことができる。コンピュータシステム７１はふつうＩ／Ｏインタフェース７４を介してディスプレイ７５ならびにマウスおよびキーボードなどのさまざまな入力装置７６に接続される。支援回路は、キャッシュ、電源、クロック回路および通信バスなどの回路を含むことができる。記憶装置７３は、ランダムアクセスメモリ（ＲＡＭ）、リードオンリメモリ（ＲＯＭ）、ディスクドライブ、テープドライブなどまたはこれらの組み合わせを含むことができる。本発明は、記憶装置７３に格納されてＣＰＵ７２によって実行されて信号源７８から信号を処理するルーチン７７として実施することが可能である。このため、コンピュータシステム７１は、本発明のルーチン７７を実行する場合に特殊用途のコンピュータシステムになる汎用コンピュータシステムなのである。
【００５３】
コンピュータシステム７１はまたオペレーティングシステムおよびマイクロ命令コードを含む。ここに記載したさまざまな処理および関数は、マイクロ命令コードの一部またはオペレーティングシステムを介して実行されるアプリケーションプログラムの一部（またはこれらの組み合わせ）とすることが可能である。さらに付加的なデータ記憶装置および印刷装置などの別のさまざまな周辺装置をこのコンピュータプラットホームに接続可能である。
【００５４】
さらに、添付の図面に示した構成システムコンポーネントおよび方法ステップのいくつかはソフトウェアで実現できるため、システムコンポーネント（または処理ステップ）間の実際の接続は、本発明がプログラムされる仕方に依存して変わり得ることを理解されたい。ここに示した本発明の教示があれば、通常の知識を有する当業者は、本発明の上記および類似の実現または構成を想到することができる。
【００５５】
本発明を有利な実施形態に基づいて詳しく説明したが、添付の特許請求の範囲に記載した本発明の精神および範囲から逸脱することなくさまざまな変更および置き換えが行えることが当業者にはわかるはずである。
【符号の説明】
【００５６】
７１  コンピュータシステム、  ７２  中央処理ユニット、  ７３  記憶装置、  ７４  入出力（Ｉ／Ｏ）インタフェース、  ７５  ディスプレイ、  ７６  入出力装置、  ７７  ルーチン、  ７８  信号源
【図面の簡単な説明】
【００５７】
【図１】（ａ）〜（ｅ）は、本発明の実施形態による顕著特徴領域の例示的な集合を示す図である。
【図２】（ａ）〜（ｄ）は、本発明の１実施形態にしたがって並進し、回転して元々のスライスに重ね合わせたＣＴボリューム（左側部分）を示す図である。
【図３】本発明の１実施形態にしたがって結合形対応付け集合およびレジストレーション変換の両方を最適化するＥＭ形アルゴリズムのフローチャートである。
【図４】（ａ）〜（ｃ）は本発明の１実施形態にしたがい、ＰＥＴ−ＣＴ，ＣＴ−ＣＴ，SPECT−ＣＴボリューム対に対してそれぞれ、ｘ，ｙおよびｚ方向においてcmで測定した全体的な距離と標準偏差とを示す表である。
【図５】本発明の１実施形態による１アルゴリズムを使用して取得され、ＰＥＴ−ＣＴ画像対，輝度アーチファクトを有するＣＴ−ＣＴ画像対，SPECT−ＣＴ画像対からそれぞれ得られ、融合されてレジストレーションされた結果画像から得られる３つのスライスを示している。
【図６】本発明の１実施形態による３Ｄレジストレーション処理のフローチャートである。
【図７】本発明の１実施形態による３Ｄレジストレーション処理を実現するコンピュータシステム例のブロック図である。
【図３】
2007054636.tif 000071
【図４】
2007054636.tif 000072
【図６】
2007054636.tif 000073
【図７】
2007054636.tif 000074
【図１】
2007054636.tif 000075
【図２】
2007054636.tif 000076
【図５】
2007054636.tif 000077