(57)【要約】【課題】事前情報を用いてマルチモードレジストレーション性能を改善すること。【解決手段】コンピュータ支援による自動画像レジストレーションの方法において、事前に位置合わせされた画像対を供給し、前記事前に位置合わせされた画像対から第１の個数の特徴サンプルを抽出し、観察された画像対から第２の個数の特徴サンプルを抽出し、前記事前に位置合わせされた画像対と前記観察された画像対とからの特徴サンプルの合併集合のユークリッド最小スパニングツリーを求め、前記最小スパニングツリーに基づいて前記観察された画像対の類似度を求め、前記類似度の勾配を推定し、前記観察された画像対のレジストレーションを出力する。【選択図】なし【特許請求の範囲】
【請求項１】
コンピュータ支援による自動画像レジストレーションの方法において、
事前に位置合わせされた画像対を供給し、
前記事前に位置合わせされた画像対から第１の個数の特徴サンプルを抽出し、
観察された画像対から第２の個数の特徴サンプルを抽出し、
前記事前に位置合わせされた画像対と前記観察された画像対とからの特徴サンプルの合併集合のユークリッド最小スパニングツリーを求め、
前記最小スパニングツリーに基づいて前記観察された画像対の類似度を求め、
前記類似度の勾配を推定し、ただし、該勾配推定値は観察された画像対のレジストレーションに使用される変換パラメータを更新するために使用されるものであり、
前記観察された画像対のレジストレーションを出力し、その際、前記観察された画像対のレジストレーションは表示されるか又は記憶媒体に記憶されるかのどちらかである、ことを特徴とするコンピュータ支援による自動画像レジストレーションの方法。
【請求項２】
前記事前に位置合わせされた画像対は、レジストレーション問題のシーケンスにおける先行する画像対であるか又はトレーニング画像対であるかのどちらかである、請求項１記載のコンピュータ支援による方法。
【請求項３】
前記ユークリッド最小スパニングツリーのパワーで重み付けされた辺の長さの勾配の和を用いて、前記類似度の勾配を推定する、請求項１記載のコンピュータ支援による方法。
【請求項４】
前記第２の個数の特徴サンプルを周期的に抽出する、請求項１記載のコンピュータ支援による方法。
【請求項５】
観察された画像対のレジストレーションを
Ｑα（Ｉ１，Ｉ２）＝Ｊα，ｗ（ｐｔ，ｐｏ）＋λ＊Ｈα（ｐｏ）
として求め、
ここで、Ｊα，ｗ（ｐｔ，ｐｏ）はα−Ｊｅｎｓｅｎ距離であり、
【数１】
2007000638.tif 000005
は項Ｈα（ｐｏ）の相対的影響を決定する重みであり、ｗ＝１−λである、請求項１記載のコンピュータ支援による方法。
【請求項６】
前記第１及び第２の個数のサンプルは同じ個数のサンプルであるか又は異なる個数のサンプルであるかのどちらかである、請求項１記載のコンピュータ支援による方法。
【請求項７】
マシンによる読み取りが可能なプログラム記憶装置であって、自動画像レジストレーションの方法ステップを実行するマシン実行可能な命令から成るプログラムを実体的に実現した形式のプログラム記憶装置において、
前記方法ステップは、
事前に位置合わせされた画像対を供給するステップ、
前記事前に位置合わせされた画像対から第１の個数の特徴サンプルを抽出するステップ、
観察された画像対から第２の個数の特徴サンプルを抽出するステップ、
前記事前に位置合わせされた画像対と前記観察された画像対とからの特徴サンプルの合併集合のユークリッド最小スパニングツリーを求めるステップ、
前記最小スパニングツリーに基づいて前記観察された画像対の類似度を求めるステップ、
前記類似度の勾配を推定するステップ、及び、
前記観察された画像対のレジストレーションを出力するステップ
を有しており、
前記勾配推定値は観察された画像対のレジストレーションに使用される変換パラメータを更新するために使用されるものである、ことを特徴とするプログラム記憶装置。
【請求項８】
前記事前に位置合わせされた画像対は、レジストレーション問題のシーケンスにおける先行する画像対であるか又はトレーニング画像対であるかのどちらかである、請求項７記載のプログラム記憶装置。
【請求項９】
前記類似度の勾配は、前記ユークリッド最小スパニングツリーのパワーで重み付けされた辺の長さの勾配の和を用いて推定される、請求項７記載のプログラム記憶装置。
【請求項１０】
前記第２の個数の特徴サンプルは周期的に抽出される、請求項７記載のプログラム記憶装置。
【請求項１１】
観察された画像対のレジストレーションは
Ｑα（Ｉ１，Ｉ２）＝Ｊα，ｗ（ｐｔ，ｐｏ）＋λ＊Ｈα（ｐｏ）
として求められ、
ここで、Ｊα，ｗ（ｐｔ，ｐｏ）はα−Ｊｅｎｓｅｎ距離であり、
【数２】
2007000638.tif 000006
は項Ｈα（ｐｏ）の相対的影響を決定する重みであり、ｗ＝１−λである、請求項７記載のプログラム記憶装置。
【請求項１２】
前記第１及び第２の個数のサンプルは同じ個数のサンプルであるか又は異なる個数のサンプルであるかのどちらかである、請求項７記載のプログラム記憶装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は画像処理に、より詳細には、事前に位置合わせされた画像対からの事前情報を用いた自動画像レジストレーションのためのシステム及び方法に関する。
【背景技術】
【０００２】
一般的な画像レジストレーションアルゴリズムは、２つの画像の間の類似性を定量化するレジストレーション測度、許容される空間変換を決定する変換空間、及び、許容される空間変換の中からレジストレーション測度を最大化（又は最小化）する変換を探す最適化スキーム、という３つの要素を持っている。
【０００３】
マルチモードレジストレーション性能の改善に事前情報を用いることに関しての問題は、Ｌｖｅｎｔｏｎ他により示唆されている。彼らは、新たな画像対に対してレジストレーション測度を決定するために、トレーニングデータと最尤法を使用して、位置合わせされた画像対の結合事前輝度分布を推定することを提案している。その後、Ｃｈｕｎｇ他が、事前に位置合わせされたデータの推定された結合輝度分布と新たな画像の結合輝度分布との間のＫｕｌｌｂａｃｋ−Ｌｅｉｂｌｅｒダイバージェンスによってレジストレーションの品質を決定する代替的アプローチを提案した。レジストレーションはこのＫ−Ｌダイバージェンスを最小化することにより達成される。ＬｅｖｅｎｔｏｎとＣｈｕｎｇの双方とも、事前情報の使用がレジストレーション方法のロバストネスを向上させることを実験的に示した。
【発明の開示】
【発明が解決しようとする課題】
【０００４】
本発明の課題は事前情報を用いてマルチモードレジストレーション性能を改善することである。
【課題を解決するための手段】
【０００５】
上記課題は、コンピュータ支援による自動画像レジストレーションの方法において、事前に位置合わせされた画像対を供給し、前記事前に位置合わせされた画像対から第１の個数の特徴サンプルを抽出し、観察された画像対から第２の個数の特徴サンプルを抽出し、前記事前に位置合わせされた画像対と前記観察された画像対とからの特徴サンプルの合併集合のユークリッド最小スパニングツリーを求め、前記最小スパニングツリーに基づいて前記観察された画像対の類似度を求め、前記類似度の勾配を推定し、ただし、前記勾配推定値は観察された画像対のレジストレーションに使用される変換パラメータを更新するために使用されるものであり、前記観察された画像対のレジストレーションを出力し、その際、前記観察された画像対のレジストレーションは表示されるか又は記憶媒体に記憶されるかのどちらかであるようにすることにより解決される。
【発明を実施するための最良の形態】
【０００６】
本発明の１つの実施形態によれば、コンピュータ支援による自動画像レジストレーションの方法は、事前に位置合わせされた画像対を供給し、前記事前に位置合わせされた画像対から第１の個数の特徴サンプルを抽出し、観察された画像対から第２の個数の特徴サンプルを抽出し、前記事前に位置合わせされた画像対と前記観察された画像対とからの特徴サンプルの合併集合のユークリッド最小スパニングツリー（ＥＭＳＴ）を求め、前記ＥＭＳＴに基づいて前記観察された画像対の類似度を求め、前記類似度の勾配を推定し、ただし、前記勾配推定値は観察された画像対のレジストレーションに使用される変換パラメータを更新するために使用されるものであり、前記観察された画像対のレジストレーションを出力し、その際、前記観察された画像対のレジストレーションは表示されるか又は記憶媒体に記憶されるかのどちらかである。
【０００７】
前記事前に位置合わせされた画像対は、レジストレーション問題のシーケンスにおける先行する画像対であるか又はトレーニング画像対であるかのどちらかである。
【０００８】
前記類似度の勾配は、前記ユークリッド最小スパニングツリーのパワーで重み付けされた辺の長さの勾配の和を用いて推定される。
【０００９】
前記第２の個数の特徴サンプルは周期的に抽出される。
【００１０】
観察された画像対のレジストレーションは
Ｑα（Ｉ１，Ｉ２）＝Ｊα，ｗ（ｐｔ，ｐｏ）＋λ＊Ｈα（ｐｏ）
として求められ、
ここで、Ｊα，ｗ（ｐｔ，ｐｏ）はα−Ｊｅｎｓｅｎ距離であり、
【００１１】
【数１】
2007000638.tif 000002
は項Ｈα（ｐｏ）の相対的影響を決定する重みであり、ｗ＝１−λである。
【００１２】
前記第１及び第２の個数のサンプルは同じ個数のサンプルであるか又は異なる個数のサンプルであるかのどちらかである。
【００１３】
本発明の１つの実施形態によれば、自動画像レジストレーションの方法ステップを実行するマシン実行可能な命令から成るプログラムを実体的に実現した形式のマシン読み取り可能なプログラム記憶装置が提供される。前記方法ステップは、事前に位置合わせされた画像対を供給するステップ、前記事前に位置合わせされた画像対から第１の個数の特徴サンプルを抽出するステップ、及び、観察された画像対から第２の個数の特徴サンプルを抽出するステップを有している。前記方法はさらに、前記事前に位置合わせされた画像対と前記観察された画像対とからの特徴サンプルの合併集合のＥＭＳＴを求めるステップ、前記ＥＭＳＴに基づいて前記観察された画像対の類似度を求めるステップ、前記類似度の勾配を推定するステップ、及び、前記観察された画像対のレジストレーションを出力するステップを有している。ただし、前記勾配推定値は観察された画像対のレジストレーションに使用される変換パラメータを更新するために使用されるものである。
【実施例】
【００１４】
本発明の１つの実施形態によれば、システム及び方法は、マルチモード医療データセット（例えば、２Ｄ及び３Ｄ、インター及びイントラ患者データ）のレジストレーションと融合とにおける自動画像レジストレーション（例えば、空間的位置合わせ）方法の性能を改善するために、事前に位置合わせされた画像対からの情報を使用する。
【００１５】
本発明の１つの実施形態によれば、方法は異なる２つのエントロピー測度、すなわち、α−Ｒｅｎｙｉエントロピーとα−Ｊｅｎｓｅｎダイバージェンスを使用する。α−Ｒｅｎｙｉエントロピー測度は、位置合わせされる画像対のみに基づいて位置合わせの品質を盲目的に定量化することを目的としている。α−Ｊｅｎｓｅｎダイバージェンスは、正しく位置合わせされた画像対のトレーニング集合からの事前情報を取り入れるために使用される。混合測度を導入することにより、これら２つのアプローチの特質が本方法に取り入れられる。
【００１６】
本発明の１つの実施形態によれば、ＥＭＳＴに基づいた混合測度の推定値を使用して、位置合わせの品質が定量化される。ＥＭＳＴは平面内の点の集合の最小スパニングツリーであり、点の各対の間の辺の重みは点の間の距離である。ＥＭＳＴに基づいた混合測度の推定値は、様々な変換のクラスに、例えば、剛体変換に適用することができ、より豊かな変換のクラスへの、例えば、自由変形、非線形パラメトリック変換モデル（薄板、Ｂスプライン）などへの拡張の可能性を示す。
【００１７】
本発明の１つの実施形態によれば、ＥＭＳＴに基づいた画像レジストレーション法において、レジストレーション測度が以下のようにして求められる。特徴対の値のリスト（スカラー、例えば、ピクセル輝度値、又は、ベクトル、例えば、ウェーブレット係数、勾配）が、位置合わせすべき画像対から抽出される。特徴リストのＥＭＳＴにおけるパワーで重み付けされた辺の長さの和が、図１Ａ−Ｃに示されているように、レジストレーション測度を決定する。図１Ａ及び１Ｂは位置合わせされるべき画像であり、ＥＭＳＴ（図１Ｃ）は結合画像からの輝度サンプルにわたって決定される。このＥＭＳＴの辺の（パワーで重み付けした）全長が位置合わせの品質の測度である。各特徴は同一の分布をする確率変数が独立して実現されたものであると仮定すれば、この方法は、数学的には、画像対のα−Ｒｅｎｙｉエントロピー（Ｈα（・））を用いることと等価である。ここで、α−Ｒｅｎｙｉエントロピーは２つの画像の間の統計的依存関係の測度であると言えることに注意されたい。レジストレーション測度は位置合わせすべき画像対からの情報のみを使用し、事前に位置合わせされた画像対からのあり得る事前情報のいずれも無視するので、ここでは盲目的なレジストレーション測度と呼ばれる。
【００１８】
本発明の１つの実施形態によれば、ＥＭＳＴレジストレーション測度において事前情報が具体化される。マルチモード画像対（Ｉ１ｋ，Ｉ２ｋ），ｋ＝０，１，２，…のシーケンスを位置合わせする状況を考えてみる。時点ｋ−１においては、画像対（Ｉ１ｊ，Ｉ２ｊ），ｊ＝０，…，ｋ−１は位置合わせされている。２つの異なる画像対の空間的位置合わせが独立していれば、事前の空間的位置合わせはＩ１ｋとＩ２ｋの位置合わせに関する直接の情報を伝えない。しかし、モダリティー関係はシーケンスに沿って不変であると仮定されているので、事前位置合わせから得られたモダリティー関係に関する情報が（Ｉ１ｋ，Ｉ２ｋ）のレジストレーションにおいて役立つ可能性がある。
【００１９】
本発明の１つの実施形態によれば、事前情報の使用がＥＭＳＴに基づいた画像レジストレーションの枠組みに取り入れられ、レジストレーション測度は事前に位置合わせされた画像対からの情報を組み込んでいる。事前に位置合わせされた画像対はレジストレーション問題のシーケンスの中の先行する画像対であってもよいし、トレーニング事例の集合からのものであってもよい。いずれの場合にも、事前情報の品質が改善されるので、同じモダリティーの間での新たなレジストレーションの正確度を改善するために、事前情報は重み付けして使用される。
【００２０】
本発明の１つの実施形態によれば、観察された分布と事前分布との間のずれに基づいて位置合わせの品質を測るために、α−Ｊｅｎｓｅｎ距離が使用される。α−Ｊｅｎｓｅｎ距離はＲｅｎｙｉエントロピーを使用して定義される確率分布間の距離測度である。α∈（０，１）とｗ∈（０，１）を固定すると、ｐｘからｐｙまでのα−Ｊｅｎｓｅｎ距離は
Ｊα，ｗ（ｐｘ，ｐｙ）＝Ｈα（ｗｐｘ＋（１−ｗ）ｐｙ）
−［ｗＨα（ｐｘ)＋（１−ｗ）Ｈα（ｐｙ）］
である。Ｈαは凹であるから、ｐｘ≠ｐｙならば、Ｊα，ｗ（ｐｘ，ｐｙ）＞０、ｐｘ＝ｐｙならば、Ｊα，ｗ（ｐｘ，ｐｙ）＝０である。ここで、ｐｘ又はｐｙはゼロに等しくない。Ｊｅｎｓｅｎ距離に基づいたレジストレーション測度の性能は事前情報の質と量に依存することになる。他方で、Ｒｅｎｙｉエントロピー測度はこの潜在的に有用な情報を無視する。本発明の１つの実施形態によれば、これら２つの測度が以下のようにして組み合わされて、混合レジストレーション測度が得られる。
【００２１】
【数２】
2007000638.tif 000003
はそれぞれトレーニング集合Ｉｔからの特徴サンプルの集合と観察された画像対Ｉｏを表し、ｐｔとｐｏは相応する分布を表すものとする。そして、混合測度を次のように定義する：
Ｑα（Ｉ１，Ｉ２）＝Ｊα，ｗ（ｐｔ，ｐｏ）＋λ＊Ｈα（ｐｏ）            （１）
ここで、
【００２２】
【数３】
2007000638.tif 000004
は、Ｒｅｎｙｉエントロピー項の相対的影響を決定する重みであり、ｗ＝１−λである。｜・｜は集合の濃度を表すことに注意されたい。
【００２３】
この測度のＥＭＳＴ推定値に基づいて、トレーニングサンプルと観察されたサンプルとを含んだ合併集合の上で求められたＥＭＳＴの辺のパワーで重み付けされた全長が、レジストレーション測度として実現される（図２Ａ及び２Ｂ参照）。
【００２４】
図２Ａ及び２Ｂはそれぞれトレーニング集合のＥＭＳＴとサンプルの合併集合のＥＭＳＴを示している。右側のＥＭＳＴの辺の（パワーで重み付けされた）全長は観察された画像の位置合わせの品質の測度である。この測度は事前に位置合わせされた画像対からの事前情報を取り込んでいる。
【００２５】
図３を参照して、以下の説明では、回転角と並進ベクトルとによってパラメトライズされた剛体変換としてパラメータ空間が定義される。３０１では、事前に位置合わせされた画像対が供給される。
【００２６】
３０２では、事前に位置合わせされた画像対から特徴のサンプルを抽出することによりトレーニング集合が形成される。特徴のタイプとしては、ピクセル輝度値、輝度勾配、及びウェーブレット係数が含まれるが、これらに限定されない。ＥＭＳＴはこの集合に基づいて決定される。このＥＭＳＴに含まれていない辺は余剰辺と呼ばれる。
【００２７】
３０３では、ある一定の個数の特徴サンプルが観察された画像対から抽出される。特徴のタイプとしては、ピクセル輝度値、輝度勾配、及びウェーブレット係数が含まれるが、これらに限定されない。正確度とスピードを改善するために、非均一サブサンプリング技術を用いてもよい。
【００２８】
３０４では、Ｋｒｕｓｋａｌアルゴリズム又は他の既存の手法を用いて、サンプルの合併集合のＥＭＳＴが求められる。これに先立って、ＥＭＳＴ法のスピードを増すために、Ｄｅｌａｕｎａｙ三角形分割を使用してもよい。さらに、余剰辺はこのＥＭＳＴを求める際に捨ててしまってよい。高次元の特徴に対して使用される択一的な方法は、（例えば、射影によって得られる）次元を下げた空間に対してＤｅｌａｕｎａｙ三角形分割を実行するものであり、この三角形分割のＥＭＳＴは原空間内で求めることができる。これにより、高次元空間内に存在しているサンプルの近似ＥＭＳＴが得られる。
【００２９】
３０５では、ＥＭＳＴに基づいた類似度が求められる。この類似度はレジストレーションの品質と変換パラメータの勾配、例えば、レジストレーションを改善するためにどのようにして変換パラメータを変更するかの指示を含んでいる。
【００３０】
抽出されたサンプルは量子化されるので、ＥＭＳＴの（パワーで重み付けされた）全長は変換パラメータに関して微分可能でない。３０６では、ＥＭＳＴのトポロジーが定正数ｋ∈Ｚ＋の回数の反復のあいだ一定であると仮定して、ＥＭＳＴメトリックの勾配の推定値、すなわち、類似度が、重み付けされた辺の長さの勾配の和として求められる。勾配推定値は観察された画像対を位置合わせするための最適化の変換パラメータを更新するために使用される。ｋ回の反復の度に、方法はブロック３０３に戻る。
【００３１】
最適化は、勾配降下法の反復の実施又はそのバリエーションであってよい。局所的な極値に閉じ込められるのを防ぐために、方法をマルチ解像度・マルチスタートの最適化と組合せてもよい。例えば、観察された画像の粗い（低解像度の）表現のレジストレーションでスタートするなど。この低解像度画像を得るために、ガウスぼけと均一サブサンプリングを使用してもよい。粗いレベルで得られたレジストレーション結果は、次のレベルでのレジストレーションの初期化に使用される。
【００３２】
３０７においていったん停止基準が満たされれば、例えば、最近Ｘ回の反復の間のパラメータ値の変化が所定の閾値を超えていなければ、方法は終了する。ここで、Ｘは所定の反復回数である。３０８では、レジストレーションされた画像対が目視のためのディスプレイ又は記憶媒体に出力される。
【００３３】
本発明はハードウェア、ソフトウェア、ファームウェア、専用プロセッサ、又はこれらの組合せの様々な形態で実施することができることを理解されたい。１つの実施形態では、プログラム記憶装置上に実体的に実現されたアプリケーションプログラムとしてソフトウェア内で本発明を実施してもよい。また、このアプリケーションプログラムを任意の適切なアーキテクチャから成るマシンにアップロードし、このマシンで実行してもよい。
【００３４】
図４を参照すると、本発明の１つの実施形態によれば、自動画像レジストレーションのためのコンピュータシステム４０１は、とりわけ、中央処理ユニット（ＣＰＵ）４０２、メモリ４０３、及び入／出力（Ｉ／Ｏ）インタフェース４０４から構成することができる。コンピュータシステム４０１は一般に、Ｉ／Ｏインタフェース４０４を介して、ディスプレイ４０５と、マウスやキーボードのような様々な入力装置４０６とに接続されている。補助回路としては、キャッシュ、電源、クロック回路、及び通信バスのような回路が挙げられる。メモリ４０３としては、ランダムアクセスメモリ（ＲＡＭ）、読み出し専用メモリ（ＲＯＭ）、ディスクドライブ、テープドライブなど、又はこれらの組合せが挙げられる。本発明をメモリ４０３に記憶されたルーチン１００７として実施し、信号源４０８からの信号を処理するＣＰＵ４０２によって実行されるようにしてもよい。コンピュータシステム４０１はそれ自体としては汎用コンピュータシステムであるが、本発明のルーチン４０７を実行する際に、特定用途向けコンピュータシステムとなる。
【００３５】
コンピュータプラットフォーム４０１はまたオペレーティングシステムとマイクロ命令コードも含んでいる。本明細書に記載された様々な処理及び関数は、マイクロ命令コードの一部であってもよいし、オペレーティングシステムを介して実行されるアプリケーションプログラムの一部（それらの組合せ）であってもよい。さらに、他の様々な周辺装置を、補助的なデータ記憶装置や印刷装置のようなコンピュータプラットフォームに接続してもよい。
【００３６】
さらに、添付図面に示されているシステムコンポーネントと方法ステップのうちのいくつかはソフトウェアで実施することもできるので、システムコンポーネント（又は処理ステップ）間の実際の接続は、本発明がどのようにプログラムされるかに応じて異なることがあるということを理解されたい。本明細書の開示内容が与えられれば、関連技術の当業者はこれらのインプリメンテーション又はコンフィギュレーション、及び、類似したインプリメンテーション又はコンフィギュレーションを想到することができるであろう。
【００３７】
自動画像レジストレーションのためのシステムと方法の実施形態を説明したが、ここで、上記の開示内容に照らして当業者により改善及び変形が為されうることを指摘しておきたい。したがって、変更は本発明の範囲及び趣旨の範囲内にある本発明の実施形態において為されるということが理解されるべきである。
【図面の簡単な説明】
【００３８】
【図１Ａ】Ｔ１強調ＭＲＩ画像を示す。
【図１Ｂ】Ｔ２強調ＭＲＩ画像を示す。
【図１Ｃ】本発明の１つの実施形態による図１Ａ及び１Ｂの結合画像からの輝度サンプルにわたって求められたＥＭＳＴを示す。
【図２Ａ】本発明の１つの実施形態によるサンプルのトレーニング集合のＥＭＳＴを示す。
【図２Ｂ】本発明の１つの実施形態によるサンプルの合併集合のＥＭＳＴを示す。
【図３】本発明による方法のフローチャートを示す。
【図４】本発明の１つの実施形態によるシステムを示す。
【符号の説明】
【００３９】
４０１  コンピュータシステム
４０４  Ｉ／Ｏインタフェース
４０７  ルーチン
【図３】
2007000638.tif 000007
【図４】
2007000638.tif 000008
【図１Ａ】
2007000638.tif 000009
【図１Ｂ】
2007000638.tif 000010
【図１Ｃ】
2007000638.tif 000011
【図２Ａ】
2007000638.tif 000012
【図２Ｂ】
2007000638.tif 000013