JP4917733B2 - 尤度最大化を利用した画像位置合わせシステム及び方法 - Google Patents

Description

【０００１】
［発明の分野］
本発明は画像処理システム及び方法に関し、特に、２枚以上の画像を１枚の合成画像に結合する画像登録システムに関する。本発明は医療用画像化の分野に特に応用されるが、本発明は、複数の画像が相互に関係し、１枚の合成画像に結合されるような、その他タイプの画像化システムにも応用できることは容易に理解されるであろう。
［発明の背景］
ここで「画像」という用語は、特に断らない限り３次元（３Ｄ）の立体的画像を指す。医学分野においては、多様な画像化方式によって立体画像を得ている。例えば、磁気共鳴像（ＭＲＩ）、Ｘ線コンピュータ断層撮影法（ＣＴ）、陽電子射出断層撮影法（ＰＥＴ）及び単光子射出断層撮影法（ＳＰＥＣＴ）などの核画像法、超音波などである。このように得られた立体画像は、ボクセル（voxel）値の配列として、通常、コンピュータのメモリなどにデジタル的に保存される。各ボクセル（voxel）は３D空間（例えばX、Y、Z座標）の位置と結びついており、色値が、典型的にはグレイスケール強度値が割り当てられている。
画像融合、あるいは複数の関連した画像のデータを集積した合成画像を形成するためのそれら画像の結合を病院環境で行うことが所望されている。多くの場合、合成した画像を見ることで、画像を別々に見ていたのでは得られない識見を診断者は得ることができる。異なる画像化方式は相補的なので、複数方式の画像融合がしばしば役に立つ。例えば、コンピュータ断層撮影法（CT）と磁気共鳴（MR）像からは基本的に解剖学的な構造上の情報が得られるのに対し、単光子射出断層撮影法（SPECT）及び陽電子射出断層撮影法（PET）からは機能上かつ代謝上の情報が得られる。機能上または代謝上の画像と、構造上または解剖学上の画像の組み合わせによって、機能上の画像の位置が特定でき、診断の精度が向上する。例えば、腫瘍学の分野では、機能上の画像における局部限定が正確に位置づけできれば、臨床医が病変の進行や治療の効果を評価することができる。こういった診断上の研究は、外科的な、あるいは目標細胞を取り囲む健康な細胞への影響を最小限に抑えなければならない放射線療法的な治療計画にも使用される。同じ方式で得られた画像を合成したい場合もある。例えば、MR血管造影像、コントラストを強調したMR像、機能的MRI像（ｆMRI）などの複数のMRスキャンの結果を解剖学的MR像などに合成することが所望される。
【０００２】
複数の画像から意味のある合成画像を形成するには、画像が適当に位置合わせされることが重要である。画像位置合わせによって画像があいまいさ無く結合できるよう、空間的な位置合わせをする必要がある。多数の画像位置合わせ技術が知られている。
【０００３】
ある画像位置合わせ技術では、その画像に表れている被写体の構造に詳しい専門家が、位置合わせする画像それぞれに目印をつける。２枚の画像は、あらかじめ分かっている目印間の関係に基づいて位置合わせされる。この方法による画像位置合わせでは、位置合わせ精度が選択された目印の数と位置によって制限されてしまう。目印の数が少な過ぎると位置合わせが不正確になる。目印の数が多くても必ずしも位置合わせが正確となるわけではなく、位置合わせに要する計算を複雑にするだけである。また、手作業の部分で時間がかかる。さらに、すべての画像で適当な構造上の目印を見出せるわけでもない。
【０００４】
最近、２種類の異なった画像化方式が１つの画像化装置に結合された。この様に画像位置合わせにおいてハードウェアを結合する手法は、コスト及び論理的な理由から、画像位置合わせの問題への最適な解答ではない。多くの場合、ハードウェアによる位置合わせは現実的ではなく、不可能でさえあり、ソフトウェアによる位置合わせ技術に依存せざるを得ない。例えば、ハードウェア手法は、例えば、時間的な治療効果を監視するとか、被験者間あるいは地域的な比較のために、異なった時間または異なった被験者から得られた画像を位置合わせすることには使えない。ソフトウェア位置合わせは、たとえハードウェアに基づく手法が採用されたとしても必要となるだろう。例えば、PETやSPECTの透過及び射出スキャンなど、同一の機械で順番にスキャンしたものの動き補正や、既定の治療計画に関して患者の位置決めをするには、ソフトウェア位置合わせが必要となろう。
【０００５】
近年、全ボリュームベースのアルゴリズムが、データの滅失に依存せず、セグメンテーション不要で、ユーザの関与がほとんど、または完全に不要であるという理由で、広く使われるようになった。さらに重要なことは、これらアルゴリズムが完全に自動化でき、位置合わせ結果を数量的に評価できることである。エントロピー・ベースのアルゴリズム、特に相互情報アプローチは、全ボリュームベースのアルゴリズム中、最も卓越したものである。これらアルゴリズムのほとんどは、２つの方式で得られた画像データを関係付ける目標関数を最適化する。エントロピーあるいは相互情報テクニックは一方向では制限されているが、別方向では制限されていない。例えば、相互情報は下限が０だが、上限は実装によって決まり、例えば、ビン（bin）の数で決まる。同様に、エントロピーは下限が０だが、上限は実装によって決まる。結果として、画像がうまく位置合わせできる程度を知ることは困難である。
【０００６】
従って、本発明は、上記その他の問題を解決する新規で改良された画像処理システム及び方法を提供することを意図するものである。
［発明の概要］
第１の態様では、本発明は、グレイスケールボクセル（voxel）値の３次元配列から成る第１と第２の立体画像を位置合わせする方法を提供する。変異確率は、前記第１の画像のボクセル（voxel）値と前記第２の画像の空間的に対応するボクセル（voxel）値から成るボクセル（voxel）値の複数の整列したペアに対して画定される。変異確率は位置合わせされた画像の統計から得るか、現在のデータの集まりから計算される。各変異確率は一方の画像のボクセル（voxel）値が他方の画像の空間的に対応するボクセル（voxel）値に対応する尤度に関係し、画像の選択された幾何学的関係に基づく。第１の画像に対して第２の画像の幾何学的関係を画定する第１の変換が選択され、変異確率を使って整列されたボクセル（voxel）値のペアの所定の集まりについて、尤度の計量が計算される。尤度の計量は、第２の画像が与えられた時に第１の画像が得られる確率と、第１の画像が与えられた時に第２の画像が得られる確率を表す。第１の画像に対して第２の画像の幾何学的関係を画定する異なる変換が選択され、尤度の最適な計量を与える最適な変換が計算されるまで、処理が繰り返される。
【０００７】
本発明の別の態様によれば、第１の立体画像と第２の立体画像を位置合わせする画像処理システムは、位置合わせプロセッサ及び位置合わせされる複数の立体画像表示を記憶する付属メモリが設けられ、前記最適な変換を表すパラメータを記憶するための前記位置合わせプロセッサに付随したメモリが設けられ、前記第１と第２の画像の合成画像表示を形成するための表示システムが設けられる。特に、前記位置合わせプロセッサは、前記ボクセル（voxel）値の複数の整列したペアの変異確率を画像の選択された幾何学的関係に基づいて決定し、前記ボクセル（voxel）値の整列した各ペアは前記第１の画像のボクセル（voxel）値と前記第２の空間的に対応したボクセル（voxel）値から成る。変異確率は前記第１の画像のボクセル（voxel）値が前記第２の画像の空間的に対応したボクセル（voxel）値に対応する尤度に関係する。前記位置合わせプロセッサは前記第１の画像に対する前記第２の画像の幾何学的関係を画定する第１の変換選択し、前記変異確率を使って整列したボクセル（voxel）値の所定の集まりについて前記尤度の計量を計算し、前記尤度の計量は前記第２の画像が与えられた時に第１の画像を得る前記確率及びその逆を表す。前記位置合わせプロセッサは異なる変換を選択し、前記尤度の最適な計量を与える最適な変換が計算されるまで処理を繰り返す。
【０００８】
本発明の利点はデータを減少せず、セグメンテーションやユーザによる操作を必要としない点である。
【０００９】
本発明の利点は確率解釈に基づく位置合わせ方法を与えることで、それは理解しやすい。
【００１０】
本発明の別の利点は位置合わせが対称的であることで、即ち第１の画像を第２の画像に位置合わせされたり、または第２の画像が第１の画像に位置合わせされたりするのではなく、２枚の画像がお互いに位置合わせされるのである。
【００１１】
本発明の別の利点は、本発明がセグメンテーションを容易に取り入れられ、それにより確率的関係が立体データの一部を使って評価可能であり、その立体データの一部は例えば、立体の空間的セグメンテーションあるいはボクセル（voxel）値に基づくセグメンテーションである。立体の一部の重要性を強調することによって、場合によってはよりよい位置合わせ結果が得られることもある。こういった柔軟性は全ボリュームをベースにした従来技術では得られなかったものである。
【００１２】
本発明の更なる利点と利益は、当業者が以下の好ましい実施形態の詳細な説明を読んで理解すれば明らかと成ろう。
［発明の実施の形態］
一般的に２枚の位置合わせされた画像のボクセル（voxel）値は確率的に関係しており、そのボクセル（voxel）値が同じ方式で得られたものか別の方式で得られたものであるかにはよらない。この確率的関係は実験的な観察に基づいて得られるし、立体データに基づいて首尾一貫して評価することができる。この確率的関係に基づいて、ある立体画像が与えられた時に、別の立体画像がある尤度を計算できる。この尤度が最大値をとる時２枚の画像を位置合わせすべきと考える。２枚の画像が位置合わせされると、正規化された最大の尤度は確率解釈が可能で、従来技術の相互情報法やエントロピー法よりずっと分かりやすい。計算された尤度の下限は０であり、上限は１となる。
【００１３】
図１を参照するに、本発明による画像処理システム１００は、立体画像データを取得し記憶するための第１の画像ソース１１０と第２の画像ソース１１２を含む。第１と第２の画像ソース１１０と１１２は、例えばMRスキャナ、X線CTスキャナ、核カメラ（PETやSPECTスキャナ等）、超音波スキャナなどの診断用画像スキャナであればよい。第１と第２の画像ソース１１０と１１２は、同じ画像方式でも異なる画像方式でもよく、異なるスキャナハードウェアで得られたものでも同一のスキャナハードウェアで得られたものであってもよい。例えば、第１と第２の画像ソース１１０と１１２は、複数の画像モードを含む１個の装置でもよい。また、第１と第２の画像ソース１１０と１１２は、異なるタイミングで複数の画像を取得する１個の装置でもよい。
【００１４】
ある実施形態では、第１と第２の画像ソース１１０と１１２は、位置合わせする画像を生成するのに適切なものとして画像診断に関する当業者に知られているセンサー、データ取得回路、画像再構成回路を含む。しかし、別の意図された実施形態においては、画像ソース１１０と１１２のいずれかまたは両方が、前もって取得された画像であり、例えば、電子的メモリまたはコンピュータで読取可能な記憶媒体に保存された画像表示であり、例えば、コンピュータメモリ、ランダムアクセスメモリ（RAM）、ディスク装置、テープ装置、その他の磁気的または光学的媒体、コンピュータネットワーク上の記憶場所などに保存される。それゆえ、本発明の画像処理システムは位置合わせする画像の一方または両方を取得したスキャナハードウェアに直接インターフェイスされるが、必ずしも直接インターフェイスされている必要はない。画像処理装置１００は画像ソース１１０と１１２からの画像データを記憶する画像メモリ１２０をさらに含む。位置合わせプロセッサ１３０は画像メモリ１２０から２枚の立体画像を読み込み、尤度最大化アルゴリズムを使って画像を位置合わせし、その２枚の画像に関する位置合わせ変換行列を作る。
【００１５】
ここで「尤度」とは、簡潔に言えば確率関係を指し、換言すると、画像ｆが与えられた時に画像ｇを得て、画像ｇが与えられた時に画像ｆを得る確率である。この尤度は複数の位置合わせまたは変換パラメータから計算可能で、尤度が最大になるように位置合わせパラメータを決定するため繰返し最適化される。従来の最適化技術を使用でき、例えば、Press et al., Numerical Recipes in C: The art of Scientific Computing, (2^nd ed.), Cambridge: Cambridge Univ. Press, 1999, Chapter 10に説明されている最適化技術を使用が使用される。最適化された位置合わせ変換行列はメモリ１５０に記憶される。
【００１６】
ある実施形態では、知識ベース実装が使われ、オプションのメモリ１４０が設けられ、画像のボクセル（voxel）ペア間の変異確率を記憶する。変異確率統計は以前の経験に基づくもので、例えば、２枚の画像の以前の位置合わせに基づくもので、尤度の計算に使用できる。２枚の画像の相対的な方向と位置を変更でき、ボクセル（voxel）ペアに基づいて尤度を計算できる。尤度を繰返し最適化することにより、例えば従来の最適化技術を使って最適化することにより、２枚の画像の最適な位置合わせを決定することができる。
【００１７】
別の実施形態では、自己矛盾のない実装が施され、変異確率についてあらかじめ知識が無くともよい。このアプローチでは、現在の位置合わせが変異確率を評価するのに使われ、この評価により尤度が計算される。２枚の画像の相対的な方向と位置を変更でき、グレイスケールペアに基づいて、新しく変異確率を計算し、尤度を計算することができる。従来の最適化技術を使って尤度を繰返し最適化して、最適な位置合わせを決定することができる。
【００１８】
随意にセグメンテーションを制限することができ、例えばユーザが入力できる。セグメンテーション制限によって選択したボクセル（voxel）について尤度を計算できる。ボクセル（voxel）は、例えば、画像の１個以上の部分ボリュームに限定できるし、ボクセル（voxel）強度値の一定の範囲に限定することもできる。その他のボクセル（voxel）の性質によるセグメンテーションも考えられる。
【００１９】
位置合わせ変換行列は、表示あるいはビデオプロセッサ１６０で使用され、画像メモリ１２０から読み出した２枚の画像を調整し、コンピュータその他の人間が読み取れる表示装置１６２上に位置合わせ変換行列により規定された合成画像としてそれらを表示する。標準的なデータ処理とプログラミングの技術を使って、画像とそれに付随する行列、及び以前の変異確率統計やセグメンテーション制限を、インデックス、ポインタ等とともに記憶する。
【００２０】
変換行列を記憶する代わりに、あるいはそれに加えて、一度最適な位置合わせを決定したら、２枚の画像から成る合成画像を記憶する。しかし、一般的には変換行列を記憶する方が好まれる。一方の画像を、位置合わせパラメータに基づいて、他方の画像の空間にフォーマットし直すこともできる。
【００２１】
画像合成以外の目的を持った画像位置合わせも考えられる。例えば、この教示による画像位置合わせは、複数の部分的に重複した画像に施すことができ、それらからより大きな１枚のボリューム画像を作成することができる。
【００２２】
ある実施形態では、位置合わせプロセッサ１３０と表示プロセッサ１６０は、従来の表示装置と組み合わせられた従来のコンピュータのソフトウェアとして実装される。
【００２３】
図２を参照するに、位置合わせプロセッサ１３０のモジュール図を示している。図１に示したように、位置合わせされる第１と第２の画像１１０と１１２が位置合わせプロセッサ１３０に入力される。知識ベース実装の場合、変異確率統計１４０が入力される。セグメンテーション制限１４２も、もしあれば入力される。位置合わせパラメータ１５０、の集まりが位置合わせプロセッサ１３０から出力されるが、位置合わせパラメータの集まりとは、即ち、変換行列の要素であり、一方の画像に適用した時、固定された座標軸系に対してその画像を変換し、２枚の画像を調整するためのものである。
【００２４】
位置合わせプロセッサ１３０はフレームメモリ１３２と１３４、尤度計算機１３６、最適化プロセッサ１３８を含む。尤度計算機１３６と最適化プロセッサ１３８は、集積処理システムのモジュールでよく、あるいは、複数のプロセッサに分散し並行処理の利益を享受できるものでもよい。動作においては、位置合わせプロセッサ１３０はフレームメモリ１３２と１３４の２枚の画像を読み取り、現在の位置合わせに関する尤度値を計算する。最適化プロセッサ１３８はその後、フレームメモリ１３４内で画像を変換し、尤度計算機１３６が尤度を再計算する。尤度計算のプロセスについては、図３に関連して以下に詳述する。尤度計算機１３６で計算される尤度地が最適になるまでステップが繰り返される。変換パラメータをメモリ１５０にアウトプットするか、ユーザに表示する。
【００２５】
もしセグメンテーション制限が入力されていれば、注目したボクセル（voxel）のみ、即ち、１以上の特定された部分ボリューム内のボクセル（voxel）、または特定範囲の値を持ったボクセル（voxel）が、尤度計算に使用され、特定された部分ボリューム外、あるいは特定範囲制限外のボクセル（voxel）は無視される。
【００２６】
図３に２枚の画像を組み合わせる典型的なプロセス３００を図示する。プロセス３００は、統合されたスキャナ／画像処理システムの一部分としてソフトウェアに実装することができ、あるいは、パーソナルコンピュータ、ワークステーション、その他の情報処理システムなどのスタンドアロンコンピュータに実装される、分離した画像処理システムでもよい。２枚の画像を組み合わせることを引例として説明しているが、そのプロセスは２枚以上の画像を組み合わせることにも容易に適用可能であり、例えば、この教示によるプロセスをすべての画像を位置合わせするまで順々に繰り返すこともできる。
【００２７】
位置合わせする２枚の画像を取得する（ステップ３０４と３０８）。画像は同一形式でもよいし、別形式でもよい。尤度の値を次のように計算する。２枚の画像f(x, y, z)とg(x, y, z)はグレイスケール(u₁, u₂, ..., u_n)と(v₁, v₂, ..., v_n)を持つと仮定する。
【００２８】
これらグレイスケールの頻度はそれぞれ(p₁, p₂, ..., p_n)と(q₁, q₂, ..., q_n)である。すると、
【００２９】
【数１】

単一形式の場合は雑音や画像化された対象の変化によって、また複数形式の場合は異なる形式が持つ内在的な特性によって、画像fのグレイスケール値u_iは画像gのグレイスケール値v_j、ただしj=1, 2, ..., mと対応する。異なる形式の画像中のボクセル（voxel）値の間には具体的な関係は仮定しないし、形式の画像内容にはいかなる制限もない。むしろ、それらの関係は純粋に統計的に記述され、知識に基づくか自己矛盾がない方法でデータだけから予測される。
【００３０】
ボクセル（voxel）u_iのv_jに対する変異確率をw_ijとする。すると、
【００３１】
【数２】

この時、pとqには次の関係が成り立つ。
【００３２】
【数３】

また、
【００３３】
【数４】

が成り立つ。
【００３４】
画像fと変異確率w_ijが与えられた時、画像gを得る条件付尤度は、
【００３５】
【数５】

ここで、積は重なり合ったボリューム全体で計算してもよいし、その部分ボリュームで計算してもよく、その場合選択された(u_i, v_j)のペアのみ、または選択された領域／ボリュームの(u_i, v_j)ペアのみが使われる。w_ij ＜1だから、この積はゼロに非常に近い値となる。計算を簡単にするため、これの対数値を使った方がよい。対数尤度は、
【００３６】
【数６】

対数尤度はある条件を満たす時最大値を取り、その条件の下では尤度も最大値を取る。上に定義した対数尤度はグレイ値ペアの数に依存する。この依存性を避けるため、正規化された対数尤度を使うほうが好ましく、
【００３７】
【数７】

ここにNはグレイ値ペアの数である。Nは重なり合ったボリューム中のすべてのペアが対数尤度の計算に使われるとすると、重なり合ったボリューム中のボクセル（voxel）数の合計である。
【００３８】
gからfへのボクセル（voxel）変異を考えると、
【外１】

を得る。この２個の対数尤度を足し合わせると（即ち、尤度を掛け合わせると）、
【００３９】
【数８】

最適な画像位置合わせを見つけるため、複数の方向について尤度を繰返し計算する最適化アルゴリズムを使用し、上に定義した対数尤度を最大化する変換を見つける。
【００４０】
特筆すべきことは、対数尤度を２つ足し合わせると、２枚の画像に関して画像位置合わせが対称的になることで、即ち、２枚の画像の順序は問題とならない。さらに、対数尤度を足し合わせると、目標機能の取り込み範囲が増大する。もし一方が特定の位置で局所的に最大値を取るなら、そして、もう一方が同じ位置で局所的に最大値を取らないならば、最適化プロセスがその位置で終了することはない。
【００４１】
もし２枚の立体画像が位置合わせされていれば、正規化された尤度は確率的な解釈をすることができる。正規化された尤度は第１の画像中のボクセル（voxel）が第２の画像中のボクセル（voxel）に対応し、第２の画像中のボクセル（voxel）が第１の画像中のボクセル（voxel）に対応する平均的な確率である。尤度は本質的に確率なので、その範囲は０から１の範囲であり、対数尤度の上限値は０である。最適化された尤度を見ると、位置合わせ品質の数量的な評価が可能である。画像が自分自身に位置合わせされれば、正規化された尤度は１となる。さすれば、形式内の位置合わせでは形式間の位置合わせの場合と比べて正規化された尤度が高く（１に近く）なる。
【００４２】
上述の通り、尤度の計算には、ある画像のボクセル（voxel）値の他の画像のボクセル（voxel）値に対する変異確率の知識が必要である。この値は２通りの方法で導き出すことができる。第１の実現は知識ベースの、即ち、変異確率に関する既知の知識に基づく。第２の実現は自己矛盾のないアプローチであり、既知の知識を前提としない。ステップ３１２において、尤度計算に使える画像に関する変異確率統計が既知か判断する。もし既知なら、処理はステップ３１６（知識ベースアプローチ）に進む。もし既知でなければ、処理は３４０（自己矛盾のないアプローチ）。
【００４３】
ステップ３１６において、以前の経験に基づき、グレイスケールペア(u_i, v_j), i=1, 2, ..., n; j=1, 2, ..., m間の変異確率の統計、または確率密度関数を有している。注目しているボクセル（voxel）ペアはスキャンされ、グレイスケールペア(u_i, v_j)及び既知の変異確率に基づいて、各ボクセル（voxel）ペアの尤度がステップ３２０で計算される。ステップ３２４において、上に説明したように、対数尤度が計算される。対数尤度は複数の異なる位置合わせに対して計算され、その計算は、対数尤度が最大化され、よって２枚の画像の最適な位置合わせが見つかるまで続く。特に、最適化アルゴリズムがステップ３２８で使われ、対数尤度が最大値になったか判断する。もし尤度がまだ最大化されていなければ、新しい変換が選択され（ステップ３３２）、２枚の画像の相対的な方向と位置が、最適化アルゴリズムによって変更され（ステップ３３６）、処理はステップ３２０に戻る。処理は対数尤度が最大化されたと判断されるまで続く（ステップ３２８）。一旦対数尤度が最大化され、最適な位置合わせが見つかると、尤度を最大とする位置合わせ変換の変換パラメータが出力される(ステップ３６４)。
【００４４】
もし変異確率が既知でなければ（ステップ３１２）、処理はステップ３４０に進み、本発明の尤度最大化プロセスは自己矛盾がない方法で実現される。「自己矛盾がない」という言葉を使うのは、現在の位置合わせによって生じた重なり合ったボリュームを使用して、尤度を計算するのに使った変異確率を評価するからである。最適化した尤度で最適化した位置合わせがなされる。
【００４５】
２枚の画像の現在の位置合わせの下では、変異確率は次のように評価することができる（ステップ３４０）。重なり合ったボリューム全体にわたってグレイ値ペア(u_i, v_j)をチェックする。Ni個のペア(u_i, v_j), j=1, 2, ..., mがある場合、(u_i, v_j)となるN_ij個のペアの中で、見積もり変異確率w_ijは次のように見積もられる、
【００４６】
【数９】

再度、重なり合ったボリューム全体、あるいはその部分ボリューム、にわたって尤度を計算する。対数尤度が重なり合うボリューム全体にわたって計算される場合、この見積もったw_ijを正規化された対数尤度に代入すると、
【００４７】
【数１０】

ここで、p_i=N_i/N及びp_iw_ij=h_ijであり、h_ijは同時確率の見積もり値である。これは否定された条件付エントロピーH(f|g)であることに注意せよ。
【００４８】
同様に、
【００４９】
【数１１】

であり、結局、
【００５０】
【数１２】

動作としては、注目しているボクセル（voxel）ペアをスキャンし、グレイスケールペア(u_i, v_j)と変異確率に基づき、ステップ３４４で各ボクセル（voxel）ペアの尤度を計算する。ステップ３４８で、対数尤度を上述のように計算する。もし重なり合うボクセル（voxel）全体にわたって尤度が計算された場合は、同時限界ヒストグラムを一くくりとする事で、前述のLLの式のように、計算を簡単化できる。対数尤度を繰返し複数の異なる位置合わせに対して計算し、新たに見積もり変異確率を新しい位置合わせに対して計算する。各位置合わせに対して、最適化アルゴリズムは現在の位置合わせが対数尤度を最大化するか判断する（ステップ３５２）。もし尤度がまだ最大とならない場合は、新しい変換を選択し（ステップ３５６）、２枚の画像の相対的な方向と位置が最適化アルゴリズムによって変更され（ステップ３６０）、処理はステップ３４０に戻る。処理は対数尤度が最大化したと判断されるまで（ステップ３５２）繰り返される。一旦対数尤度が最大化され、最適な位置合わせが見つかったら、尤度を最大化する位置合わせ変換の変換パラメータを出力する（ステップ３６４）。
【００５１】
知識に基づく実装及び自己矛盾のない実装両方の場合に、計算された最大化された尤度そのものを、例えば表示モニターに出力する（ステップ３６４）。計算値は上限値１の数値なので（例えば、画像を自分自身と位置合わせする場合、あるいは１対１に対応した２枚の画像の場合）、その値は画像が相互にいかによく位置合わせされたかを示す目安となる。この点、結果の値がとくに制限されないエントロピー・ベースや相互情報ベースとは異なる。
【００５２】
再述するに、知識に基づく実装であっても自己矛盾のない実装であっても、重なり合うボリューム全体にわたって尤度を計算する必要はない。ある実施形態では、所定の点、あるいはセグメントできる重なり合うボリュームの部分について尤度を計算することができる。ボクセル（voxel）値を強調することもでき、例えば、値が特定の範囲、例えば５０％以上、２５％以上、１５％以上等にあるボクセル（voxel）を強調することができる。特定のセグメントを強調することによって、位置合わせ精度を上げることができる。
経験的試験
本発明による尤度最大化画像位置合わせ方法の精度と信頼性を次の試験で実証した。特に、自己矛盾のない実装が使われ、重なり合ったボリューム中のグレイ値ペアが変異確率を評価し対数尤度を計算する。
【００５３】
これらの例では剛体変換を使ったが、より一般的な変換を考えることもでき、例えば、非線形変換、アフィン変換、写像変換等を考えることもできる。剛体変換は、位置合わせパラメータは６次元ベクトル(θ_x,θ_y,θ_ｚ, t_x, t_y, t_z)であり、θ_x、θ_y、θ_zはそれぞれｘ、ｙ、ｚ軸の周りの度で表した回転角であり、t_x、t_y、t_zはそれぞれｘ、ｙ、ｚ軸方向のｍｍで表した移動オフセットである。各回転について、等質の座標系について４ｘ４行列が対応している。さらに回転を施すと、それは行列の積で表せる。行列の積は交換可能ではないので、回転の順序が重要である。仮定によって、回転はｘ、ｙ、ｚ軸の周りで、この順番に施されるとする。さらに仮定により、平行移動の前に回転が施されるものとする。
【００５４】
自己矛盾のない実装では、重なり合ったボリューム中でグレイ値ペアの限界的同時分布を評価する。画像ｆの最大ボクセル（voxel）値をまず見つける。画像ｆのボクセル（voxel）値をn_f個の離散値に分割する。同様に画像ｇのボクセル（voxel）値はn_g個の離散値に分割する。n_fとn_gは違う個数でよい。重なり合うボリュームにおいて、画像ｆとｇのボクセル値とボクセルペアのヒストグラムがボクセル値とボクセル値ペアを組みにする（binning）事で計算される。ｆのヒストグラムのビン（bin）数はn_f、ｇにヒストグラムのビン（bin）数はn_g、同時ヒストグラムのビン（bin）数はn_fｘn_gである。以下に説明するように、ビン（bin）の数は固定することも、動的に変更することもできる。正規化したヒストグラムは同時分布であると同時に限界分布を示す。
【００５５】
変換を施した後には、一方のボリュームのグリッド点は変換された空間における他方のグリッド点と普通一致しない。ボクセル値とボクセル値ペアを組にする（binning）する前に、補間を施して変換した空間のグリッドにおけるボクセル値を得なければならない。補間法にはたくさんの異なった方法があるが、例えば、最近接法、３線法、３線部分ボリューム法等がある。１ボクセルまでの並進移動には影響を与えないので、最近接法はボクセル以下の精度が必要な時には充分ではない。簡便なので、３線法が好まれる。
【００５６】
変換の下、複数次元方向セット最適化を使用し、符号を反転した対数尤度を最小化（尤度を最大化）する。方向行列はユニタリー行列に初期化される。ベクトルは上述の通り、(θ_x, θ_y, θ_z, t_x, t_y, t_z)である。これら位置合わせパラメータの順序を変えることもでき、より速く最適化できるかも知れない。パラメータの順序の最適化は試みていないが、順序は画像内容に依存するかも知れないし、すべてを試みるのは現実的ではないからである（６！＝７２０通りの組み合わせがあり、この一部だけを試みるにしても現実的ではない）。さらに、最適化には探索が進むにつれて意図していた６個の方向とは必ずしも対応しないさらに６個の独立な方向が使われるかも知れない。
【００５７】
真に全体最適な値を見つけるため、シュミレーテッドアニーリングを利用することもできる。シュミレーテッドアニーリングは２次元画像の位置合わせには成功裏に適用できた。確率的な方法で速度も遅いので、３次元画像の位置合わせが限界であろう。実際には、複数解またはサブサンプリングアプローチが役に立つだろう。強いアルゴリズムで最適化のスピードを上げ、取得範囲を拡大する。
【００５８】
この実装においては複数解最適化が好ましい。画像は最も粗い画像として８ｘ８ｘ８の画像に折りたたまれる。３方向すべてで画像解像度が最高になるまで、続く画像の解像度を２倍ずつしていく。ボリュームが５１２ｘ５１２ｘ５１２等と大きいと、計算が大変で、位置合わせに時間がかかる。ある場合には、最高解像度の画像は使わない。解像度を低くするため、サンプリングボリューム内のボクセル値を平均化する。サブサンプリングアプローチよりは少し遅いが、実際には位置合わせ結果はよくなる。
【００５９】
同時２次元ヒストグラムを評価する時、グレイ値は異なったビン（bin）でペアとする。同時分布は正規化された２次元ヒストグラムで評価されるので、統計的観点から、サンプルサイズは大きいほうがよい。複数解最適化という戦略においては、低解像度画像が使われるが、グレイ値ペアの数は少ない。統計的バイアスが大きいと思われる。
【００６０】
画像サイズを８ｘ８ｘ８としよう。その時、最大で５１２個のグレイ値ペアがある（すべてのボクセルが重なり合う時）。８ビットのグレイデータでは、ビン数は最大２５６となる。ビンが２５６個あっては、平均して各ビンに最大２個のペアしかなくよくない。同時確率の統計エラーによって結果が悪くなる。この状況では、ビン数は少ないほうがよい。
【００６１】
ビン数を小さな値に固定しておけば、精細な画像の場合に、各ビンに充分なグレイ値ペアがある。感度を犠牲にして同時確率の評価をよくすることもできる。このパラドックスはビンの数を固定するのは理想的ではなく、ビン数を適合させる方がよい、即ち解像度に応じてビンの数を変化させる方がよいことを暗示している。
【００６２】
テストでは、ビンの数を適合させた。ビン数が多いと位置合わせ処理が遅くなるので、最大ビン数は１２８とした。ビンの適合数を解像度の関数として表１に示したが、ここで解像度は３方向におけるピクセルの最大数で与えられている。
【００６３】
【表１】

本発明による尤度最大化画像位置合わせ法を次のように評価した。SPECTとMRのデータセットをテストボリュームとして使用した。SPECT画像とMR画像に様々な位置ずらしをかけた後に、ボクセル値にノイズを入れた場合と入れない場合で、自己位置合わせを行った。そして、SPECTのボリュームをMRのボリュームに位置合わせした。位置合わせの結果を評価した。
【００６４】
画像データはスライスから成る。ｘ軸は右から左に水平に、ｙ軸は前から後ろへ水平に、ｚ軸は下から上に垂直に伸びている。テクネチウム−９９ｍヘキサメチル−プロピレネアミン−オキシム（HMPAO TC-99m）をSPECT画像の取得用の試薬として使った。画像は６４ｘ６４ｘ２４ボクセルのサイズで、ボクセルサイズは７．１２ｘ７．１２ｘ７．１２ｍｍ^３であった。最小ボクセルグレイ値は０で、最大ボクセル値は５４２５であった。平均ボクセル値は１６３．８であった。
【００６５】
MR画像（T1矢じり状の）は２５６ｘ２５６ｘ１２８ボクセルで、ボクセルサイズは１．０ｘ１．０ｘ１．５ｍｍ^３であった。最小ボクセル値は０で最大ボクセル値は５０４、ボクセル値の平均は５７．８であった。
【００６６】
まず、様々な位置ずれから始めて、画像をそれ自身に位置合わせした。ノイズに対する耐性を評価するため、一方の画像を異なったレベルのノイズでよごし、ランダムに作った位置ずれの同一セットで位置合わせを行った。ランダムに位置ずれした画像ペアを作るため、画像をまずｘ軸の周りにある角度回転し、ｙ軸の周りに別の角度回転し、ｚ軸の周りにさらに別の角度回転した。これらの角度は一定の範囲で一様な分布をしている。回転した画像を新しい位置に動かした。ｘ、ｙ、ｚ方向のオフセットはある範囲で一様に分布している。
【００６７】
画像にホワイトノイズを加え、画像をよごした。ノイズはガウス分布で、平均値は０、標準偏差は様々に変化させた。
【００６８】
これらの素性の知れた実験によって、位置合わせ結果を検査した。熟練の観察者だからｘ軸、ｙ軸方向に２ｍｍ以上、ｚ軸方向に３ｍｍ以上の並進位置ずれを発見でき、ｚ軸の回りに２°以上、ｘ軸、ｙ軸の回りに４°以上の回転位置ずれを認め、もし位置ずれパラメータの内どれか一つでもこれらの範囲を超えると、位置合わせは失敗であると考えた。５０組のランダムに位置ずれを起こしているボリュームのペアを位置合わせした。これらの位置ずれを起こさせるために、一方の画像をｘ、ｙ、ｚ軸の回りに回転した。この回転角は−２０°から２０°の範囲で一様に分布させた。その後、画像を並進移動させた。３軸に沿った並進オフセットは−５６．９６から５６．９６ｍｍの範囲で一様に分布させた。
【００６９】
ノイズのない核画像については、本技術による方法は５０回中１０回だけ失敗した。比較のため、これら位置をずらしたボリュームを相互情報最大化法及び平均絶対差最小化法によっても位置合わせを行った。相互情報アプローチは５０回中４１回失敗し、平均差アプローチは５０回中４回の失敗であった。
【００７０】
位置合わせを成功させるため、平均回転ずれと並進ずれが、位置合わせパラメータの標準偏差とともに計算し表２に表示してある。核画像は異なったレベルのノイズと位置ずれを含む自分自身と位置合わせした。各欄において、第１の数字は平均値、第２の数字は標準偏差である。角度は度、並進はｍｍ、時間は秒で表している。３種類のアルゴリズム、平均絶対差最小化法（AD）、相互情報最大化法（MI）、尤度最大化法（LH）を評価した。
【００７１】
【表２】

位置ずれパラメータは実際の値と計算値との違いとして定義される。それらはボクセル以下の精度で得られ、バイアスと偏差が非常に少ない。本発明の尤度最大化法は相互情報アプローチより約２５％高速だった。画像を位置合わせした時、計算した尤度は正に１で、ボクセルが１対１対応していることを示している。
【００７２】
本発明による尤度最大化法の強さを評価するために、画像を付加的ノイズで汚した。同一の初期位置ずれで再度画像を位置合わせした。位置ずれパラメータの統計は、位置合わせが成功したものについて表２に表示してある。ガウスノイズの分散は５０と１００とした。平均差法は付加的ノイズがガウスノイズでないと失敗するが、これらの方法はN(0, 100)というノイズがあってもうまく動いたことは特筆すべきである。これだけ強かったのは、画像信号が強かったからだろう。平均グレイ値は１６３．８であった。ノイズがN(0, 50)の場合、解像度が６４ｘ６４ｘ２４の時、位置合わせした時の計算された尤度は０．７６２、解像度が３２ｘ３２ｘ２４の時、０．９４４であった。同様に、ノイズがN(0, 100)の場合、尤度はそれぞれ０．５９４と０．９０３であった。ボクセル範囲は非常に広いので、各ビン（bin）は比較的多くのボクセルグレイ値を持つことができた。その結果、ノイズで汚れたボクセルが、ノイズフリーのボクセルとして同じビンに入ったのかも知れない。それが、尤度が大きくなった理由である。低い解像度の画像を使った場合、ボクセルを滑らかにし、ビンの数を減らしたことが、低い解像度で高い尤度が得られることを説明している。ここのデータはノイズが高ければ尤度が低くなることを示している。
【００７３】
同じ実験をMR像についても行った。しかし、回転角度は−３０から３０度の範囲にわたって一様に分布しており、並進オフセットは−２４から２４ｍｍの範囲で一様に分布している。
【００７４】
MRボリュームは２５６ｘ２５６ｘ１２８ボクセルから成る。前に述べたように、全解像度を使っている場合は、位置合わせの速度は遅い（約１時間）。従って、最大解像度は６４ｘ６４ｘ６４とした。次のデータが示すとおり大変うまくいくことが分かった。
【００７５】
表３は、位置ずれにおいてノイズがある場合とない場合とについて、成功した位置合わせにおける位置ずれパラメータの統計を示している。各欄で、第１の数字は平均値で、第２の数字は標準偏差である。角度は度で、並進移動はｍｍで、時間は秒で表している。３つのアルゴリズム、平均絶対差最小化法（AD）、相互情報最大化法（MI）、尤度最大化法（LH）を評価した。成功率はほぼ同じオーダーだが、相互情報はすべての画像ペアを正しく位置合わせした。位置合わせ結果への解像度の影響は調べなかったが、それは失敗した位置合わせはすべて大きな並進オフセットを持っており、重なり合うボリュームがないなどが理由である。この状況はより高い解像度の最適化でも変わらない。
【００７６】
【表３】

ノイズフリーの画像を位置合わせする場合、尤度は予想通り１であった。ノイズを含む画像の位置合わせをする場合、解像度が１２８ｘ１２８ｘ１２８の時、尤度は０．０８２、解像度が６４ｘ６４ｘ６４の時、０．２４８であった。前述の例のように、解像度が下がると尤度が上がる。ここで尤度が前の例より小さいのは、平均ボクセルグレイ値がここでは低く、ノイズの影響がより大きかったからである。
【００７７】
複数形式の位置合わせでは、正しい位置合わせパラメータは分からない。多様な評価方法がその正確性を評価するために使用されているが、その中には、ファントム確認、オブザーバ評価、基準マーク等がある。ここでは、モディファイドオブザーバ評価アプローチを使用した。MR／SPECT像ペアは最初医療専門家によりソフトウェアとして手に入るインターラクティブ（マニュアル）位置合わせ法を使って位置合わせされた。位置合わせパラメータはその後基準として使われ、尤度最大化法または相互情報最小化法の位置合わせ結果が、この基準結果と比較された。これら技術の耐久性を評価するため、最初の位置合わせの時に、位置ずれをランダムに起こした。５０セットの位置ずれ（セット１）において、回転角度と基準回転角度の差は−２０から２０度の範囲にわたって一様に分布しており、並進オフセット間の差は−２８．４８から２８．４８ｍｍの範囲にわたって、一様に分布している。５０セットの第２の位置ずれ（セット２）について、角度の差は−３０から３０度の範囲で一様に分布しており、並進オフセットの差は−２８．４８から２８．４８ｍｍの範囲で一様に分布している。真の位置合わせパラメータは分からないので、比較的大きな閾値を計算値と基準値の間に設け、位置合わせが成功したか判断した。角度の閾値は回転軸に関わらず１０度で、並進の閾値は方向に関わらず１４．２４ｍｍであった。
【００７８】
最初、MR像は基準として使われ、SPECTを浮遊像とした。相互情報最大化と尤度最大化の位置ずれパラメータの統計は表４に示す通りである。
【００７９】
【表４】

相互情報最大化アプローチの成功率は本発明による尤度法よりわずかによかった。当初の位置ずれが大きくなるにつれ、両方法の成功率は低下した。データから、最初の位置ずれが大きくなるにつれ、統計はあまり変わらなく、両方法とも耐性があり信頼性があることが分かった。
【００８０】
両方の位置合わせ結果は、標準のものよりシステマティックな差異があるが、それらはすべてボクセル以下である。相互情報アプローチの偏差は尤度アプローチの偏差より小さく、その差はごく僅かであるが、相互情報位置合わせはわずかに堅実である。
【００８１】
それほど大きな差は期待できないが、SPECT像を基準として使った時、実験結果は表５に示した様に、不均衡を示した。
【００８２】
【表５】

尤度アプローチの成功率は大きく低下した。これら２つのアプローチの位置合わせパラメータの偏差も拡大した。両方のアルゴリズムにおいて、位置合わせ結果はどちらの画像を基準にしたかにはよらないはずである。この不均衡は実装によるものと思われる。
【００８３】
位置合わせにおいて、尤度は解像度が３２ｘ３２ｘ２４の時、０．０６０で、解像度が６４ｘ６４ｘ２４の時、０．０２４であった。これらの値は大きさが小さいように見えるが、尤度としては大きな値である。いかに議論するように、もし関係がないなら、尤度はそれぞれ０．０００９８と０．０００２４となるだろう。
【００８４】
本発明の自己矛盾のない実装において、重なり合うボリューム全体にわたって尤度を計算する場合、最後の式はエントロピーと同時エントロピーからの貢献によって成る。導出において、画像ｆをランダムなフィールドとして取り扱えば、uiはランダム変数となる。条件付尤度を計算する代わりに、尤度を計算できる。式は符号を反転させた同時エントロピーとなる。画像が所与なので、条件付尤度を計算するほうが好ましい。前述の通り、対数尤度は符号を反転させた条件付エントロピーであるが、その解釈はまったく異なる。もし部分ボリュームデータを使って尤度を計算するなら、あるいは尤度が既知の知識に基づいて計算されるなら、式は異なり、エントロピー、条件付エントロピー、相互情報アプローチには対応するものはない。
【００８５】
しかし、エントロピーアプローチには概念的な違いがある。エントロピーの概念は１９４０年代に提唱された。エントロピーは従来、統計物理、評価理論、スペクトル分析、画像再生、キュー理論、都市計画等で最大化される目的関数として利用されている。１９８０年代の中ごろ、エントロピー最大化は公理系に基づく原理として確立された。しかし、原理としてのエントロピーの最小化は未だ提案されていない。画像位置合わせにおいてエントロピー最小化をインフォーマルにかつ直感的に正当化することはできるが、適切な基礎付けがまだなされていない。しかし、尤度最大化は確立された原理であり理解もしやすい。尤度最大化は強固な理論的基礎付けを有し、概念的に理解しやすい。
【００８６】
尤度は確率であり、下限は０、上限は１である。尤度の計算は解像度とビン数にも依存するが、上限はあくまでも１で変わらない。２枚の画像が相関していなければ、ｆ中のボクセルのｇ中のボクセルに対する変異確率は1/n_fであり、ｇ中のボクセルのｆ中のボクセルに対する変異確率は1/n_fである。その時、尤度は1/n_fn_gである。尤度を見ることで、２枚の画像がどの程度マッチしているかがわかる。エントロピーと相互情報は絶対的な概念であるが、尤度はより容易に理解できる。
【００８７】
よって、尤度最大化による画像位置合わせを使った画像処理システム及び方法を説明した。相互情報最大化法及び絶対差最小化法を比較して、上記結果は本発明の性能は相互情報アプローチの性能と同程度であることを示している。
【００８８】
本発明をある程度詳細に説明したが、要素は当業者によって本発明の精神と範囲を逸脱することなく変更できることを認識すべきである。本発明の実施形態のひとつは１またはそれ以上のコンピュータシステムのメインメモリにある命令のセットとして実現可能である。コンピュータシステムが要求するまでその命令のセットは別のコンピュータ読取可能なメモリに記憶されていてもよく、例えば、ハードディスクや、CD-ROMまたはDVDドライブで使用する光ディスク、フロッピディスクドライブで使うフロッピディスク、フロプティカルドライブで使うフロプティカルディスク、またはパーソナルコンピュータのカードスロットで使うパーソナルコンピュータメモリカードなどのリムーバブルメモリに記憶しておいてよい。さらに、命令のセットは別のコンピュータのメモリに記憶しておき、ローカルエリアネットワークやインターネットなどのワイドエリアネットワークを通して、ユーザが欲する時に伝送するようにしてもよい。加えるに、命令はアプレットの形でネットワークを通して伝送してもよく、その場合、伝送前ではなくコンピュータに伝送後に解釈される。当業者は理解しているであろうが、命令またはアプレットのセットの物理的な記憶は、それが電気的、磁気的、化学的、物理的、光学的、あるいはホログラフとして記憶されている媒体によって、物理的に変化し、媒体がコンピュータによって読取可能な情報を格納できればよい。
【図面の簡単な説明】
本発明は多様な構成要素及び構成要素の配列から成り、多様なステップ及びステップの配列から成る。図面は好ましい実施形態を例示する目的で描かれたもので、本発明を限定する意味に解釈してはならない。
【図１】 本発明による画像取得及び処理システムのブロック図である。
【図２】 本発明による立体画像位置合わせ（registration）プログラムを実装した場合の多様なモジュールを例示するブロック図である。
【図３】 本発明による画像位置合わせ（registration）処理のフローチャートである。

Claims (20)

1. 第１の立体画像と第２の立体画像を位置合わせする方法であって、各画像はグレイスケールのボクセル値の３次元配列を含み、
   （ａ）複数の前記ボクセル値の整列したペアの変異確率を画定するステップであって、前記ボクセル値の整列したペアは前記第１の画像からのボクセル値と前記第２の画像からの空間的に対応したボクセル値とを含み、各変異確率は前記第１の画像のボクセル値が前記第２の画像の空間的に対応したボクセル値に対応する尤度及びその逆に関係し、前記画定する段階は前記第１の画像と前記第２の画像の選択された幾何学的関係に基づくステップと、
   （ｂ）前記第１の画像に対する前記第２の画像の幾何学的関係を画定する第１の変換を選択するステップと、
   （ｃ）前記変異確率を用いて整列したボクセルのペアの所定の集まりの尤度の計量であって、前記第２の画像を前提として前記第１の画像を得る確率及びその逆を表す尤度の計量を計算するステップと、
   （ｄ）前記第１の画像に対する前記第２の画像の幾何学的関係を画定する異なる変換を選択するステップと、
   （ｅ）前記第１の画像に対する前記第２の画像の幾何学的関係を画定する最適な変換であって前記尤度の計量が最適になるものが求まるまで（ｃ）と（ｄ）のステップを反復するステップとを含む方法。
2. 前記最適な変換を表すデータを記憶するステップ、及び
   前記最適な変換を用いて前記第１と第２の画像を位置合わせするステップの少なくとも一方をさらに含む、請求項１に記載の方法。
3. 前記第１と第２の画像から形成される合成画像を表示するステップをさらに含む、請求項２に記載の方法。
4. 前記最適な変換は剛体変換、アフィン変換、写像変換、非線形変換から選択される、請求項１乃至３いずれか一項に記載の方法。
5. 前記変異確率を画定するステップはあらかじめ計算され記憶された変異確率値にアクセスするステップを含む、請求項１乃至４いずれか一項に記載の方法。
6. 前記変異確率を画定するステップは、前記第１の画像に対する前記第２の画像の幾何学的関係を画定する現在選択されている変換に基づいて変異確率を推定するステップを含む、請求項１乃至５いずれか一項に記載の方法。
7. 前記複数の整列したペアは前記第１と第２の画像の重なった部分にあるすべての整列したペアを含む、請求項１乃至６いずれか一項に記載の方法。
8. 前記複数のペアは前記第１と第２の画像の重なった部分の１以上の選択された部分ボリュームに限定される、請求項１乃至７いずれか一項に記載の方法。
9. 前記複数のペアは、前記第１と第２の画像の重なった部分にある整列したペアで、所定範囲のボクセル値を持つペアに制限される、請求項１乃至８いずれか一項に記載の方法。
10. 前記尤度の計量は対称的である、請求項１乃至９いずれか一項に記載の方法。
11. 前記尤度の計量を計算するステップは前記変異確率の対数を合計するステップを含む、請求項１乃至１０いずれか一項に記載の方法。
12. 前記尤度の計量を計算するステップは、
    前記変異確率の対数の合計として、前記第２の画像に対して前記第１の画像を得る前記確率を表す第１の対数尤度値を、前記所定の整列したボクセルペアの集まりについて計算するステップと、
    前記変異確率の対数の合計として、前記第１の画像に対して前記第２の画像を得る前記確率を表す第２の対数尤度値を、前記所定の整列したボクセルペアの集まりについて計算するステップと、
    前記第１と第２の対数尤度値を足し合わせるステップとを含む、請求項１乃至１１いずれか一項に記載の方法。
13. 前記第１と第２の対数尤度値を足し合わせるステップは所定範囲の合計尤度を与え、さらに、位置合わせ品質が上がるにつれ前記合計尤度値が前記範囲内で上がる、請求項１２に記載の方法。
14. 前記尤度の正規化された計量を作るため、前記変異確率を使って、整列したボクセルペアの所定の集まりの尤度の前記計算された計量を正規化するステップであって、前記尤度の前記正規化された計量は、前記第２の画像に対する前記第１の画像を得る確率、及び前記第１の画像に対する前記第２の画像を得る確率を表し、前記正規化された計量は所定の有限の範囲に制限されているステップと、
    前記正規化された計量を出力するステップであって、前記位置合わせの前記品質が向上すると前記正規化された計量が前記有限の範囲の一方の端に近づき、前記正規化された計量は前記位置合わせの前記品質をユーザに示すように働く、請求項１乃至１３いずれか一項に記載の方法。
15. 第１の立体画像と第２の立体画像を位置合わせする画像処理システムであって、前記立体画像はボクセル値の３次元配列を含み、
    位置合わせプロセッサ及び位置合わせされる複数の立体画像表示を記憶する付属メモリを有し、
    前記位置合わせプロセッサは、前記ボクセル値の複数の整列したペアの変異確率を決定し、前記ボクセル値の整列したペアは前記第１の画像のボクセル値と前記第２の画像の空間的に対応したボクセル値を含み、各変異確率は前記第１の画像のボクセル値が前記第２の画像の空間的に対応したボクセル値に対応する尤度及びその逆に関係し、前記決定は前記第１と第２の画像の選択された幾何学的関係に基づいて決定され、
    前記位置合わせプロセッサは、前記第１と第２の画像の複数の幾何学的関係の尤度の計量を計算し、前記尤度の計量は前記変異確率を使って整列したボクセルペアの所定の集まりについて計算され、前記尤度の計量は前記第２の画像に対する第１の画像を得る前記確率及びその逆を表し、
    前記位置合わせプロセッサは、前記第１と第２の画像の幾何学的関係を画定する最適な変換を見つけるため前記尤度の計量を最適化し、
    前記画像処理システムは、さらに
    前記最適な変換を表すパラメータを記憶する、前記位置合わせプロセッサに結合したメモリと、
    前記第１と第２の画像の合成画像表示を形成する表示システムとを有する、画像処理システム。
16. 診断画像スキャナを更に有する、請求項１５に記載の画像処理システム。
17. 前記診断画像スキャナはMRスキャナ、X線CTスキャナ、PETスキャナ、SPECTスキャナ、超音波スキャナ、またはこれらの組み合わせである、請求項１６に記載の画像処理システム。
18. 前記変異確率の決定は記憶された計算された変異確率値にアクセスすることを含む、請求項１５乃至１７いずれか一項に記載の画像処理システム。
19. 前記変異確率の画定は、前記第１の画像に対して前記第２の画像の幾何学的関係を画定する現在選択されている変換に基づき変異確率を推定することを含む、請求項１５乃至１８いずれか一項に記載の画像処理システム。
20. 前記尤度の前記計量は対称的である、請求項１５乃至１９いずれか一項に記載の画像処理システム。

Images