対象物の移動による画像アーチファクトの減少を伴う生物医学画像の登録のための方法と、コンピュータプログラムプロダクト
登録されるべき第１と第２画像において、一定の数の特徴又はランドマークが自動的に定義され、追跡されて、第１と第２画像間の光学フローベクトルが定義されることを特徴とする、対象物の移動による画像アーチファクトの減少を伴う生物医学画像の登録のための方法が開示される。登録は、逆光学フローを第２画像に加えることで行われる。定義された各画素に対する隣接領域における平均信号強度を定義し、平均信号強度を所定の閾値と比較することによって、自動特徴選択ステップが実行される。もし前記近隣領域の平均信号強度が、所定の閾値より高ければ、その画素は特徴として定義され、追跡すべき特徴のリストに加えられる。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、
ａ）画素の2次元又は３次元アレイによって形成された画像からなる、同一対象物の少なくとも第１と第２デジタル画像又は横断面画像セットを提供するステップと、
ｂ）前記第１画像又は画像セットの中で、一定数のランドマーク又は所謂特徴を、一定数の画素を選択することによって定義し、追跡すべき特徴のリストを生成するステップと、
ｃ）前記第１又は第２画像又は画像セットの中から特徴として選択された画素を有する画像の光学的フローベクターが定義されることによって、特徴として選択された各画素の位置が追跡されるステップと、
ｄ）逆光学フローを第二画像又は画像セットに加えることによって、第一と第二画像又は画像セットが登録されるステップとからなる、画像の登録のための方法に関する。
【背景技術】
【０００２】
部分的に時系列のこの種の画像の登録方法は、例えばＵＳ６５５３１５２によって知られている。
【０００３】
ＵＳ６５５３１５２文献は、人間オペレーターによる視覚的ランドマークの識別と両画像の対応する画素のマーキングによる２画像の画像登録法を提供している。
【０００４】
この方法は、ランドマークの選択をオペレーターの技能に完全に委ねており、画像に著しく明確な認識可能な構造が存在しない場合に、実行することは実用的に非常に困難である。例えば磁気共鳴画像（MRI）又は超音波又は放射線画像といった診断用画像を考慮した場合、画像の構造的特徴によるランドマークの視覚的識別は、極めて困難であり、大きなエラーが生じる可能性がある。
【０００５】
様々な異なるアルゴリズムが輪郭を識別しようとしてきた。例えば、Fischer, B . ; Modersitzki, J. Curvature Based Registration with Applications to MRMammography; Springer-Verlag Berlin Heidelberg: ICCS 2002, LNCS 2331, 2002; pp 202- 206に開示されている。又は、2回以上測定された組織／領域の放射線画像内のランドマークや、異なる方法で調査されたものとして、例えばShi J, Tomasi C . Good Features to Track . 1994 IEEE Conference on Computer Vision and Pattern Recognition (CVPR' 94) 1994; 593- 600 and in Tommasini T, Fusiello A, Trucco E, Roberto V. Making Good Features Track Better Proc . IEEE Int . Conf . on Computer Vision and Pattern Recognition (CVPR ' 98) 1998; 178- 183 and in Ruiter, N. V. ; Muller, T . O . ; Stotzka, R. ; Kaiser, W. A. Elastic registration of x- ray mammograms and three-dimensional magnetic resonance imaging data . J Digit Imaging 2001, 14 , 52-55.が開示されている。
【０００６】
異なる時間に同一対象物の2枚の平面画像を登録し、対象物の時間に対する移動量を計算することを、2次元登録と呼ぶ。同様の処理を3次元の画像セット、例えばMRIやコンピューター断層撮影（CT）スキャンの横断面画像などに用いるアルゴリズムを、３次元と呼ぶ。ランドマークや対象物の移動は、これら画像セットの中でいろんな方向に生じる。
【０００７】
比較すべき画像の異なる様式の数に依存して、アルゴリズムは単一モード登録と、マルチモード登録アルゴリズムとに分かれる。
【０００８】
２枚の連続したMRI画像を比較することは、単一モード登録である。
【０００９】
通常、単一モード登録アルゴリズムは、ほとんどの場合2次元で働く両画像内の一定のランドマーク又は複雑な形状を識別しようとする。第二ステップにおいて、両画像内の対応するランドマークの位置が比較され、それらの移動ベクトルが計算される。この情報は、第２画像のすべての画素を、新しい位置に移動させ、移動アーチファクトを除去するために用いられる。
【００１０】
「厳格な登録」は、単一ユニットとしての2次元又は３次元画像を特定の方向に移動する。「非厳格な登録」は、一定領域の単一又は複数の画素を異なる方向に移動する。
【００１１】
例えば非常にしなやかな胸の組織に関し、非厳格な登録が移動アーチファクトを除去するために必要である。
【００１２】
骨のように形成された構造の無い特別な一貫性のために、直接隣接した胸の部分は、異なる方向に移動可能である。従って、いかなる登録アルゴリズムも他の実在するアプローチも、胸の組織上のランドマークを識別し、広げることは明らかである。
【００１３】
一つのアルゴリズムグループは、有限数の発見すべきランドマークを定義する。これらは、それらの価値に従って記憶される。胸のMR画像に関し、MRIの脂肪と空気間の高いコントラストのおかげで、これらのアプローチは胸の外境界又は胸部構造部で多くの数のランドマークを見つけることがある（Huwer, S . ; Rahmel, J. ; Wangenheim, A. v . Data-Driven Registration for Lacal Deformations . Pattern Recognition Letters 1996, 11, 951-957）。
【００１４】
圧縮パネルを用いるか、又は少なくともうつ伏せの位置で横になることで両胸を固定することで、通常は圧縮装置に直接接触する肌のすべての部分を良好に固定することができる（図１）。対照的に、両胸の内部領域は、心臓の鼓動と呼吸によって僅かに動き続ける。脂肪と柔組織間の外側皮膚と空気間と比較して非常に低いコントラストにより、これらのアルゴリズムにおいて胸の中心にある有効なランドマークの数は、低く留まっている。
【００１５】
ランドマークを選択することにおける一般的な問題点は、非常に少ないランドマークの選択により、不十分な又は不正確な登録が行われることである。しかし、追加のランドマークを選択することは、正確な登録を保証するために必要ではない。なぜなら、人間の観察者又は何らかのプログラムが、例えば胸の中央の低いコントラストや複雑な組織タイプによって生じる再割り当てされた低い確実性のランドマークを含めるように強いられるからである。ランドマーク数の上昇は、常に登録計算の複雑さを著しく増加させる。
【００１６】
所謂非厳格な登録に特に適した特徴選択と特徴追跡アルゴリズムは、所謂ルーカス＆カナダアルゴリズムと、その特有のピラミッド型実行であり、詳細には以下に開示されている（Lucas BD, Kanade T . An Iterative Image Registration Technique with an Application to Stereo Vision , IJCAI81; 674-679 1981 and in Bouguet J- Y . Pyramidal Implementation of the Lucas Kanade Feature Tracker, TechReport as Part of Open Computer Vision Library Documentation, Intel Corporation）。これらは付録１として下に同封する。ルーカス＆カナダ技術と、そのピラミッド型特徴追跡実行は、特に画像中の確かなランドマークの自動識別や、異なる時間での二つの画像間のランドマークの追跡に適している。選択された特徴の置換や移動は、移動量又は所謂光学的フローベクトルとして測定される。この技術は、例えばロボット工学の画像認識などに広く応用されており、デジタル画像処理の当業者には一般的な知識の一部であり、画像認識はいくつかの基本コースで教えられている。
【００１７】
しかし、デジタル診断画像、より具体的には胸画像へのルーカス＆カナダ技術の応用の実用的な用途を考えると、この技術は、興味の外側まで含めた非常に多数の特徴を識別し、例えば胸の中の空気といったノイズの影響を受ける。
【００１８】
以下の例では、上記問題のより詳細な観点を与える。
【発明の開示】
【発明が解決しようとする課題】
【００１９】
コントラストを強調した（CE）磁気共鳴乳房撮影法（MRM）は、侵略する腫瘍に対し素晴らしい感度を有する有用で重要な調査法である。それは、胸の柔組織の密度と独立な高い負の予測値を示す。最近のマルチセンター研究では、感度と特異性は、画像解析技術に従って異なることが示されている。1.5Tスキャナーを用いた、９６、９３又は８６％の感度で３０、５０又は７０％の特異性のものがそれぞれ計算されている。従って。MRMの主な欠点は、誤った好結果の診断が高い確率で生じることである。
【００２０】
通常MRイメージングは、1.0又は1.5テスラ撮像装置でメーカーの二重ブレストコイルを用いて（米国では、しばしば単一ブレストコイルが用いられる）、患者をうつ伏せにして行われる。プロトコルは、少なくともあらゆる向きで捕捉時間３０〜１８０秒でのダイナミックコントラストで強調されT1で重み付けされたシーケンスを有する。少なくとも一つの測定が、ガドペンテト酸ジメグルミンの大丸薬注入又は他の常磁性のMRコントラスト薬剤の投入前に行われ、幾つかの測定がその後に行われる。
【００２１】
胸の組織の領域への実際のコントラスト薬剤の取り込みを計算するために、コントラスト薬剤投与前の各画素の信号強度は、コントラスト薬剤投与後の信号強度から引き算されなければならない。MR画像のスライス厚さで生じる部分的な体積効果と同様に、呼吸又は心臓の鼓動によって生じる両胸の最小限の移動によって、異なる時間に得られた同じ画像位置の画素は、正確に胸の組織の同じ厚さを示さない。従って引き算された画像は、絶対的に黒にはならない（腫瘍部を除く）。少しの移動の効果は、胸の外側境界を良く説明することができるであろう。脂肪組織の高い信号強度と周囲の空気のおよそゼロの信号強度とによって、非常に小さな移動が、脂肪を空気の以前の画像位置に置く。結果として、引き算の画像は、画像化された胸の少なくとも部分周りに取り込まれた強いコントラスト薬剤を現わす厚い白い線を示す（図２で、胸の白い外側境界がサインによってマークされている）。
【００２２】
デカルト体系によって定義された3次元又は体積画像を得ることを考慮すると、ｘ、ｙ方向へのあらゆる移動に加えて、（胸に向かう）ｚ方向へのいくらかの移動が常に見られる。これは、一連の調査中の患者のくつろぎによって生じる。
【００２３】
アーチファクトを減らすために、胸は通常ある種の圧縮法（メーカーに依存する）によって胸コイルの内側に固定される。しかし、この事実にかかわらず、最小の移動は常に発生する。従って、胸を固定しているにもかかわらず、画像には僅かな移動アーチファクトがｘ、ｙ、ｚ方向に見られる。移動アーチファクトは明るい画素として引き算画像の中で視覚化される。移動アーチファクトがなければ、引き算画像は腫瘍組織に占領された領域を除き、絶対的に黒である。
【００２４】
一定数のランドマークが胸の外側の周囲の空気のノイズの大きな信号の中から検出される。考慮された各特徴が追跡されるべきであり、不必要な計算上の負荷が発生するという事実の他に、ノイズはランダムな効果であって、ノイズに含まれる全ての特徴又はランドマークは、第２の、又は後に続く画像の対応するランドマークを有さない。しかし事実として、ノイズ内で見出されたランドマークに対応するアルゴリズムがサーチする第２画像やそれに続く画像内のランドマークは、ともかく第２画像や画像セット内に一定量のランドマークを再割り当てする。これが、画像の間違った登録に繋がり、結果に影響を与える。
【課題を解決するための手段】
【００２５】
本発明の目的は、画像、特に3次元生物医学的診断用画像を登録する方法に存する。この方法は、3次元での移動追跡を可能とし、デジタル画像内の特に興味のある領域の外側のノイズによるアーチファクトを除去し、隣接する組織の信号強度の高低による差異と独立に胸にかかるすべてのランドマークを広げることによって、既知の登録方法の欠点を効果的に克服するものである。
【００２６】
本発明によれば、画像を登録する方法は以下のステップからなる。
ａ）画素の2次元又は３次元アレイによって形成された画像からなる、解剖学上の部分の組織又は組織帯をＭＲＩ、超音波、又は放射線で走査して得られた第１デジタル画像又は横断面画像セットを提供するステップと、
第二の時間に同じ解剖学上の部分の組織又は組織帯をＭＲＩ、超音波、又は放射線で走査して得られた少なくとも第二デジタル画像又は横断面画像セットを提供するステップと、
ｂ）前記第１画像又は画像セットの中で、一定数のランドマーク又は所謂特徴を、一定数の画素を選択することによって定義し、追跡すべき特徴のリストを生成するステップと、
ｃ）前記第１又は第２画像又は画像セットの中から特徴として選択された画素を有する画像の光学的フローベクターが定義されることによって、特徴として選択された各画素の位置が追跡されるステップと、
ｄ）逆光学フローを第二画像又は画像セットに加えることによって、第一と第二画像又は画像セットが登録されるステップと、
ｅ）ステップｄ）での登録後に、第２画像の第２画像データから第１画像の画像データを引き算するステップと、
ｆ）第２画像の第２画像データから第１画像の画像データを引き算して得られた画像データを表示するステップ。
【００２７】
ステップｂ）の特徴選択は、以下により詳細に記載する。
Ｂ１）第一画像又は第一断面画像セットの各画素周りの隣接する画素領域が定義され、各画素領域が制限された数の画素を有し、単一画像の場合、２次元領域が選択され、断面画像セットの場合、３次元領域が選択されるステップと、
Ｂ２）各画素に対して、前記画素領域の全画素の平均信号強度が求められるステップと、
Ｂ３）平均信号強度閾値が求められるステップと、
Ｂ４）ステップ２で求められた各画素領域の平均信号強度と、前記平均信号強度閾値とが比較されるステップと、
Ｂ５）前記領域の平均信号強度が、ステップ４の閾値より大きい場合、画素は追跡されるべきランドマークとして定義され、追跡されるべきランドマークリストに追加されるステップとから成る自動特徴選択ステップ。
【００２８】
平均信号強度閾値は、調整可能であり、変更可能である。かかる閾値を決めるために、経験的な基準を適用することができる。例えばテスト画像又はサンプル画像に本方法を実行し、異なる閾値に対して得られた結果を比較する。
【００２９】
本発明によれば、画像の追跡すべき特徴を定義できるすべての種類の画素を選択する方法を適用可能であるが、最良の結果は既知のルーカス＆カナダ法による自動特徴選択法によって得られ、この方法の詳細は、Lucas BD, Kanade T . An Iterative Image Registration Technique with an Application to Stereo Vision, IJCAI81; 674-679 1981に記載されている。
【００３０】
所謂光学的フローという移動又は置換ベクトルによって、第１画像から第２画像への特徴の移動を追跡することは、既知の方法である。またこの場合、本発明をランドマークや特徴を利用する幾つかの追跡方法に適用することができる。特に、もし診断画像撮影中に制限された胸の小さな移動に関わるものである場合、そしてルーカス＆カナダ特徴選択法を使用する場合、本発明が、好ましくは所謂「ルーカス＆カナダ特徴追跡のピラミッド型実行」（付録１参照）を使用する。この「実行」は、当業者には既知であり、幾つかの大学の教育課程で教えられている。このルーカス＆カナダ特徴追跡のピラミッド型実行の特定の、詳細な内容は、Bouguet J-Y . Pyramidal Implementation of the Lucas Kanade Feature Tracker, Tech Report as Part of Open Computer Vision Library Documentation, Intel Corporation and in Shi J, Tomasi C . Good Features to Track . 1994 IEEE Conference on Computer Vision and Pattern Recognition (CVPR' 94) 1994; 593- 600 and in Tommasini T, Fusiello A, Trucco E, Roberto V. Making Good Features Track Better . Proc . IEEE Int . Conf . on Computer Vision and Pattern Recognition (CVPR ' 98) 1998; 178-183.に開示されている。
【００３１】
ひとたび全ての特徴に対して光学的フローベクトルが定義されると、逆フローを第２画像に加えることで、第２画像の各特徴とその周辺画素とを第１画像の特徴位置に移動することによって、２つの画像の登録が実行される。
【００３２】
既に述べたように、上記方法は、磁気共鳴イメージング（MRI）又はコンピュータ断層撮影法（CT）の同じ又は異なる種類のスキャナーを用いた、２次元又は３次元画像を登録するために用いることができる。
【００３３】
上述の登録方法は、MRI又は超音波画像といったコントラスト媒体強調診断イメージング、特にコントラスト強調磁気共鳴乳房X線撮影の方法を組み合わせることでも提供される。この方法において、胸の組織の領域への実際のコントラスト薬剤取り込みを計算するために、コントラスト薬剤塗布前の各画素の信号強度が、コントラスト薬剤塗布後の信号強度から引き算されなければならない。このように、画像化された組織においてコントラスト薬剤が存在しない時間における第１画像と、コントラスト薬剤が存在する第２の又は後の時間の第２画像とを登録することと、画像化された組織内を一面に覆うことが、もし登録が実行されなかった場合、患者の移動によって画像データの引き算に生じるアーチファクトを、高い度合いで抑圧するために役立つ。主となる考えは、２枚の画像間を引き算することで、コントラスト薬剤の溜まった組織に対応した画素が高い信号強度を有することである。コントラスト薬剤の存在しないその他の部分は、低い信号強度の画素として、特に暗い灰色又は黒色で表示される。
【発明を実施するための最良の形態】
【００３４】
本発明の一つの観点は、ルーカス＆カナダアルゴリズムの２次元から３次元への拡張であるが、以下のグラフィック例は、理解を容易にするために２次元画像を用いている。３次元画像に対する本方法の一般化は、数学的な形式論を用いて概要によって与えられる。
【００３５】
３次元と３次元データを現す断面画像セットの自動特徴選択作業は、ルーカス＆カナダの２次元特徴検出アルゴリズム実体的な拡張であって、以下に詳細に説明する。
【００３６】
時間を少し変えて記録した同じ患者の２枚の画像を、IとJと定義する。
【００３７】
最初の画像Iにおいて、２枚の画像間で追跡すべき特徴であると考えられる一定数の選択された画素が識別されなければならない。
【００３８】
Iをオリジナル／最初の画像とし、Jを次の画像とし、ここで対応する特徴は見出されていなければならない。
【００３９】
第一のステップは、最初の画像の最小固有値の計算からなる。
【００４０】
これは、以下のように実行される。
第１画像を、３つの勾配量Ix、Iy、Izに変換し、（デカルト座標系のｘ、ｙ、ｚ方向で画像を描写する）、一方ｘ、ｙ、ｚは各単一画素の座標である。
勾配量は、隣接する画素に対する各画素の強度勾配から算出され、隣接する画素は一般的にx方向で（２？ｘ＋１）のサイズ、ｙ方向で（２？ｙ＋１）のサイズ、ｚ方向で（２？ｚ＋１）のサイズを有する所謂ウインドウからなる。ここで、？ｘ、？ｙ、？ｚ＝１、２、３、４、．．．、ｎ 画素である。
【００４１】
目標画素の中央の３×３×３アレイの隣接画素のサイズを考慮して、３つの勾配量は以下のように定義される。
ａ）ｘ方向の勾配量
ｂ）ｙ方向の勾配量
ｃ）ｚ方向の勾配量
各画素の次のステップとして、傾斜マトリックスＧが、上述の３つの計算された勾配量を用いて表される。
傾斜マトリックスは、以下のように生成される。
各画素の各勾配マトリックスに対し、最小固有値？minが、以下に従って傾斜マトリックスより算出される。
勾配マトリックスGの固有値は、以下のように計算される。
閾値は、固有値？ｍの最小値として定義される。
画素が追跡可能な特徴を現しているか否かを考慮するために、以下の基準が適用される。
１．もし？ｍが閾値より大きい場合、画素の実際の位置が追跡すべき特徴リストに保存される。
２．もし？ｍが固有値？ｍの全ての最小値に対する最大値のパーセンテージより小さい場合、追跡すべき特徴リストから落とされる。
３．実際の特徴の周りの定義された領域（調整可能な、例えば３×３×３）に他の特徴が存在する場合、固有値？ｍのより小さい最小値を有する特徴は追跡すべき特徴リストから落とされる。
４．もし特徴位置周りの３ｄ隣接領域（例えば？＝１）の平均信号強度が、ゼロに近い非常に低い値の（例えば１０）平均信号強度閾値より小さいか又は等しい場合、この特徴は、胸の外側の周辺の空気に位置すると見なされ、追跡すべき特徴リストから落とされる。
５．最後に、定義可能な最大数の特徴が保持される。もし追跡すべき特徴リストに、それ以上の特徴があった場合、より小さい最小値？ｍを有する特徴がリストから落とされる。
【００４２】
これらのステップは、２次元の場合を参照して図３ないし６で視覚的に説明する。この場合、画素アレイは小さな四角のアレイとして示されている。異なる画素は、第一基準を満足する固有値の最小値を有する選択された特徴の位置に対応して示される。これらの特徴と位置は、対応する画素に異なる様相を与えることによって、図３で強調表示されている。図４において、５×５の隣接領域が、その内側に円を書くことで強調表示されている。
【００４３】
図１において、患者の胸を固定する典型的な様相が図示されている。通常、MRイメージングは、患者がうつ伏せ状態で、メーカーの２重胸コイルを用いて（しばしば米国においては単一の胸コイル装置が用いられる）、1.0から1.5テスライメージャーを用いて行われる。プロトコルは、少なくともダイナミックコントラスト強調T1重み付けシーケンスをすべての方向に対して利得時間３０〜１８０秒で有する（個々のプロトコルで異なり、脂肪抑圧を付加的に使用可能である）。少なくとも一つの測定が、0.1mmol/Kgから0.2 mmol/kg体重のガドペンテト酸ジメグルミンの大丸薬投与又は他の常磁性MRコントラスト薬剤投与前に、幾つかの測定が投与後に行われる。
【００４４】
少なくとも第一MRI画像は、コントラスト媒体が肘前の静脈に投与され、胸組織に達する前の、時間t0の時に得られる。少なくとももう一つの画像が、コントラスト媒体投与後の時間t1に得られる。通常、時間を遅らせた一連の画像が得られる（例えば６回以上連続して同じ対象に対しコントラスト媒体投与後に得る）。
【００４５】
現在の画像を単純化するために、断面図の3次元画像から取り出した２枚の画像は、一定の(スライス)厚さを有する。従って、2次元画素に関して論じる代わりに、3次元画素に関してより詳細に論じる。各画素は、対応する信号強度と相関関係のある所定のグレーレベルを有する。図３ないし１０の例において、グリッドは画素アレイを表現している。グリッドの各要素は、画像を形成する画素である。
【００４６】
図３において、例として４種類の画素がランドマークとして特徴づけられており、対応する勾配マトリックスの最小固有値に従って保存される。102と示された画素は、最大固有値を有する画素であり、202と示された画素は、大きな固有値を有する画素であり、302は小さい固有値を有する画素を、402は最小固有値を有する画素を示す。
【００４７】
図４において、各画素に対して直径５画素のEucliadian距離内の隣接領域が特徴として選択された中央画素の周りに定義される。この隣接画素の選択は、円として図４に示されている。
【００４８】
図４において、このイベントの第一の場合が502'で示された厚い円で強調表示されている。ここで、特徴として選択された２画素は、画素のEucliadian距離として上で定義されたと同じ隣接領域内にある。これらの2画素は、202'と302'で示され、上記定義に従って、特徴202'として選択された画素は、第二の画素302'より大きな固有値の勾配マトリックスを有する。従って、この第二画素302'は除かれる。
【００４９】
図５は次のステップを示し、図３で特徴として選択された302'はもはや存在しない。図５において、画素202'はその隣接領域を示す円を備えている。図４と同様に、202''の隣接領域に202''より小さな固有値を有する傾斜マトリックスの画素302''が存在する。
【００５０】
従って、この場合も特徴として選択された画素302''は追跡すべき特徴のリストから除かれる。
【００５１】
図６の画素アレイにおいて、少なくとも幾つかの特徴として選択された当初からの画素が存在する。単一の隣接領域の中では、一つ以上の特徴はもはや見出すことができない。従って追跡すべき特徴は、画素全体に広がった。
【００５２】
本発明の更なるステップによれば、画素の強度は更に調査され、所定値より低い強度を有する特徴として選択された画素もまた、追跡すべき特徴のリストから除かれる。ほぼゼロの非常に低い閾値が、周辺の空気内で見つけられる特徴の特徴づけに用いられる。このステップは図示されないが、当業者には自明である。
【００５３】
ひとたび画素が2次元または３次元画像の特徴として考慮すべきと決定すると、特徴の追跡が実行される。好ましき実施例によれば、これは所謂「ルーカス・カナデ特徴追跡のピラミッド型実行」を用いて実行される。詳細は、Jean Yves Bouguet in the article "Pyramidal Implementation of the Lucas Kanade Feature Tracker, Description of the Algorithm" published By Intel Corporation, Microprocessor Research Labsに記載されている。
【００５４】
上述の出版物において2次元例として開示されている理論を、以降3次元の場合に更に適用する。以下のステップは、上述の方法とアルゴリズムに従って個性化された追跡すべき特徴を現す各画素に対して実行されなければならない。
【００５５】
ボリュームIのポイントｕと、ボリュームＪの対応するポイントｖを定義する。第一ステップにおいて、二つのピラミッドが画像ＩとＪに対しレベルＬ＝０、．．．、ｍを現すＩLとＪLとして、また当初ボリュームを表現する基部（Ｌ＝０）を有して組み立てられなければならない。各階のボリュームは、直接隣接する領域の一定量の画素を一つの平均値と組み合わせることによってサイズが減少する。追加の階の数は、ボリュームのサイズに依存する。続くステップにおいて、所謂グローバルピラミッド型移動ベクトルｇが最も高い階で初期化される。
【００５６】
全てのレベルＬに対して、結果としての移動ベクトルｄLが計算される。
【００５７】
これは以下のステップに従って、実行される。
ａ）各ボリュームｌLにおいて、ポイントｕLの位置が以下の式で計算される。
ここでｐは実際のボリューム位置である。
【００５８】
ｂ）ボリュームのピラミッド型の表現の各レベルＬで、勾配はデカルト座標の各方向ｘ、ｙ、ｚに対して計算される。
【００５９】
デカルト座標ｘに対し、ボリューム勾配が下式に従って計算される。
【００６０】
この式は、他の二つのデカルト座標ｙとｚに対しても同様に適用され、
【００６１】
ｃ）ボリューム勾配は、次に以下の式に従って勾配マトリックスを計算するために用いられる。
【００６２】
ここで、値？は、画素に影響を与える隣接領域のエリアサイズを定義する。
【００６３】
マトリックスの要素は上述の定義に従った方法と同様に、下記式で得られる。
【００６４】
ｄ）各レベルＬに対して、相互作用ベクトル？が初期化される
【００６５】
【００６６】
以下の式に従って、ミスマッチベクトルが計算される。
【００６７】
光学的フローベクトルは、下記式によって定義される。
【００６８】
相互作用移動ベクトルは上で定義した光学的フローベクトルを用いて下のように計算できる。
【００６９】
そして、これがレベルＬに対する最終光学的フローと等しいセットである。
【００７０】
上のステップは、最終レベルＬ＝０まで各レベルで繰り返される。
【００７１】
次のレベルＬ−１に対し、グローバル光学フローは上記計算値から計算できる。
【００７２】
ｆ）最終光学的フローベクトルは、次式となる。
【００７３】
ｇ）ボリュームＩのポイントＵに対応するボリュームＪのポイントｖの座標は、下記のように計算できる。
ｖ＝ｕ＋ｄ
【００７４】
図７ないし１０は、現実の表現を明らかに単純化した上述の追跡方法のグラフィカルで模式的なアイディアを与えるものである。
【００７５】
図７ないし１０において、異なる時間Ｔ０とＴ１に撮影された二つの画像ＩとＪが、図３ないし６の例と同様に、画素の2次元アレイとして現されている。
【００７６】
上述の方法に従って追跡されるべく選択された特徴に対応する画素２０が、画像Ｉの２次元画素アレイ内に示されている。同じ位置の画像Ｊにおいて、画素２０’’は画像化された対象、この場合胸が動かない場合に画素２０と対応する。画素２０’’の隣接する特定の画素が、中央が画素２０’の７×７の画素アレイに対応してグレイの影付き四角で示されており、境界部は番号１０で示されている。
【００７７】
２０’は、画像Ｉの特徴２０の新たな位置を示している。これは、画素２０の位置にある胸の組織が時間と共に画素２０’の位置に動いたことを意味する。明らかに、画素２０’は追跡が実行されたことを示すためにのみ図示されている。実際の場合、この画素は知られず、又は追跡が不要である。
【００７８】
数学的な形式論によって既に述べられた選択された特徴に対応する画素を追跡する一般的な方法は、以下に簡単のために2次元の場合に限定して図面とともに実際的な例を記載する。しかし、当業者は追跡の一般的な数学的記載がそれでも理解でき、3次元の場合に追跡する方法を評価する。
【００７９】
追跡の初めに、画像Ｉの画素２０は、同じ位置、即ち画像Ｊの画素２０’’にあるものと仮定する。この画素２０’’は以降、被追跡画素と呼び、特徴として選択された一方画像Ｉの画素２０は特徴画素と呼ぶこととする。第一に、被追跡画素２０’’の隣接領域に対応する一定領域の逆勾配マトリックスにアルゴリズムを適用する。本例では画像Ｊの画素２０’’を中心とした７×７画素アレイ（この値は調整可能）を計算する。所謂追跡領域にこの隣接領域は対応する。ルーカス＆カナデの特徴追跡アルゴリズムのピラミッド型実行を適用することによって、3次元画像にまで適用でき、その詳細は後述するが、画像Ｊの追跡された特徴の新しい位置を見つけることができる。
【００８０】
画像Ｊにおける画素２０’’の追跡される特徴の新しい位置は、光学的フローベクトルを決定する。これは、第一の反復ステップで計算されることによって定義される。図７と図８の比較から明らかなように、追跡領域は、同じ組織領域２０’の新しい位置により近い追跡領域の中央になるように移動する。
【００８１】
追跡は、画像Ｉの特徴を現す画素２０と、画像Ｊの追跡領域中央の新たに識別された位置との間でミスマッチベクトルの計算を繰り返すことによって、さらに反復的に適用される。これが新しい光学的フローベクトルの定義をもたらし、図９に示すように画素隣接領域又は追跡領域の新たな移動をもたらす。
【００８２】
画像Ｉの特徴画素２０と選択領域の中央との間のミスマッチベクトルを決定するステップは、ミスマッチベクトルが調整可能な最小値に達するまで行われる。
【００８３】
この状態は、図１０に示される。図７ないし９の図で表現されているものに関し、画素２０’はその位置を維持しており、追跡領域が移動しているという事実を強調している。追跡領域は、最終的に正しい新しい位置２０’を中心とした位置に移動し、画像Ｉの画素２０としての胸の同じ位置を表現している。
【００８４】
上述の例は、ただ一つの特徴が画素によって表現されている場合に実行される。上述のステップは、画像内で選択された各特徴に対して実行されなければならない。
【００８５】
全ての特徴がその新しい位置まで追跡された後、画像上の全ての異なる位置にある、胸の組織を表現する領域内にある一定量の光学的フローベクトルが、存在する。他の登録アルゴリズムと比較して最大の長所は、このアルゴリズムが画像全体に光学的フローベクトルを広げることである。
【００８６】
次に、逆光学的フローベクトルを画像Ｊの各特徴とその周囲の隣接領域とに加えることで、モーフィングを実行する。これによって、新しい視覚的画像である画像Ｊ’となる。
【００８７】
図１１は、この方法の実際の診断画像での効果を説明している。左側の第一の画像セットは、コントラスト薬剤投与無しでの時間t0での胸の代表的なスライスを示す。左側の第二画像セットは、コントラスト薬剤投与後の時間t1での同じスライス部位を示す。左側の第三画像セットは、時間t1での各画素の信号強度から、時間t0での強度を引き算した、所謂引き算画像の結果を視覚化したものである。時間間隔t1−t0の間に相対的な胸の移動の無い場合は、結果画像は画素の広い部分でコントラスト薬剤強調の生じない黒で示される信号強度ゼロの部分が含まれる。時間間隔t1−t0間の移動は、組織の画像の同じ領域とは対応しない画像ＩとＪの同じ位置にアーチファクトが生成されるという効果を得る。本例において、引き算画像は小さな矢印で示されたコントラスト強調された腫瘍を含む一つの胸を示している。ほとんどの種類の胸の悪性と良性の腫瘍は、血管分布の増加によって生じる追加的なコントラスト薬剤の取り込みを示す。両方の胸の上部の追加的な白い構造体は、移動アーチファクトによるものである。両方の胸の中央に残った白い構造体は、移動との組み合わせによる通常の胸の柔組織の弱いコントラスト薬剤の取り込みの結果である。
【００８８】
図１１の右側は、登録アルゴリズムが適用された対応する画像を示している。
【００８９】
時間t0での画像は、アルゴリズムによっては変更されないため、同じままである。しかし、3次元特徴追跡後は、新たな視覚的画像が計算される。時間t1'での右側の画像は、登録適用後の対応する同じスライス位置のバーチャル画像を示す。右側の最後の画像は、時間t1'での各画素の信号強度から時間t0での信号強度を引き算した、画素対画素の引き算結果を新たな（バーチャル）引き算画像として視覚化したものである。一つの胸は、小さな矢で示すコントラスト強調腫瘍を示している。通常の胸の柔組織を現す胸の範囲内の領域では、いくつかの中間のコントラスト薬剤取り込みが見られ、移動アーチファクトの減少により少し暗く見えている。
【００９０】
図１１の例は、コントラスト薬剤投与前後の2枚の連続画像のみにより簡略化したものである。
【００９１】
図２から図１１の画像はすべて、同じMRスキャナで撮影し、登録以外の加工は行っていない。
【００９２】
実際には、連続した断面写真（例えば50枚以上）がコントラスト媒体の注入前後の通常の時間間隔で得られ、引き算画像がそれぞれの連続画像の組み合わせに対して計算される。
【００９３】
一連の画像の中で一枚以上の横断面画像を得るのに加えて、2連以上（例えば５以上）の画像の中から次々に得ることができる。これは、幾つかの引き算画像となる。もちろん、開示したアルゴリズムは、コントラスト薬剤移動の前後に関わらず、連続写真の全ての組み合わせに適応でき、その間に生じたいかなる移動も減じることができる。
【００９４】
以下の別表において、ルーカス＆カナデの追跡アルゴリズムのピラミッド型実行の詳細な理論を示す。
【００９５】
本発明に係る上述の方法は、コンピューターで実行可能なプログラムとしてコード化されたアルゴリズムとして実行される。
【００９６】
このプログラムは、例えば磁気的、光磁気的、又は光学的基板、例えば1枚以上のフレキシブルディスク、1枚以上の光磁気ディスク、1枚以上のCD-ROM又はDVD-ROM、その他代替的なものなどの携帯型又は静止型記憶装置や、携帯型又は据え置き型ハードディスクなどに保存可能である。
【００９７】
擬似コードで書かれたコード化されたプログラム例を以下に示す。
【図面の簡単な説明】
【００９８】
【図１】胸部MRIで用いられる胸固定手段の概要図である。
【図２】引き算過程と、登録ステップが実行されない場合の移動アーチファクトの発生を示す図である。左上と右上の図は、同じMRスキャナ位置で撮った２枚の胸の断面図である。第２画像は、第１画像撮影の後、約２．５分のコントラスト媒体投与後に測定した。両画像は、断面画像の３次元セットの代表部分である。
【図３】模式的な例と共に特徴選択基準が示された模式的な２次元画像である。
【図４】模式的な例と共に特徴選択基準が示された模式的な２次元画像である。
【図５】模式的な例と共に特徴選択基準が示された模式的な２次元画像である。
【図６】模式的な例と共に特徴選択基準が示された模式的な２次元画像である。
【図７】どのように特徴追跡が実行されるかを模式的に示した図である。
【図８】どのように特徴追跡が実行されるかを模式的に示した図である。
【図９】どのように特徴追跡が実行されるかを模式的に示した図である。
【図１０】どのように特徴追跡が実行されるかを模式的に示した図である。
【図１１】コントラスト媒体投与前後の異なる時間に撮影した胸の２枚の画像と、特徴追跡と登録実行後の引き算をした画像であり、左側には比較として図２の画像例を示す。
【技術分野】
【０００１】
本発明は、
ａ）画素の2次元又は３次元アレイによって形成された画像からなる、同一対象物の少なくとも第１と第２デジタル画像又は横断面画像セットを提供するステップと、
ｂ）前記第１画像又は画像セットの中で、一定数のランドマーク又は所謂特徴を、一定数の画素を選択することによって定義し、追跡すべき特徴のリストを生成するステップと、
ｃ）前記第１又は第２画像又は画像セットの中から特徴として選択された画素を有する画像の光学的フローベクターが定義されることによって、特徴として選択された各画素の位置が追跡されるステップと、
ｄ）逆光学フローを第二画像又は画像セットに加えることによって、第一と第二画像又は画像セットが登録されるステップとからなる、画像の登録のための方法に関する。
【背景技術】
【０００２】
部分的に時系列のこの種の画像の登録方法は、例えばＵＳ６５５３１５２によって知られている。
【０００３】
ＵＳ６５５３１５２文献は、人間オペレーターによる視覚的ランドマークの識別と両画像の対応する画素のマーキングによる２画像の画像登録法を提供している。
【０００４】
この方法は、ランドマークの選択をオペレーターの技能に完全に委ねており、画像に著しく明確な認識可能な構造が存在しない場合に、実行することは実用的に非常に困難である。例えば磁気共鳴画像（MRI）又は超音波又は放射線画像といった診断用画像を考慮した場合、画像の構造的特徴によるランドマークの視覚的識別は、極めて困難であり、大きなエラーが生じる可能性がある。
【０００５】
様々な異なるアルゴリズムが輪郭を識別しようとしてきた。例えば、Fischer, B . ; Modersitzki, J. Curvature Based Registration with Applications to MRMammography; Springer-Verlag Berlin Heidelberg: ICCS 2002, LNCS 2331, 2002; pp 202- 206に開示されている。又は、2回以上測定された組織／領域の放射線画像内のランドマークや、異なる方法で調査されたものとして、例えばShi J, Tomasi C . Good Features to Track . 1994 IEEE Conference on Computer Vision and Pattern Recognition (CVPR' 94) 1994; 593- 600 and in Tommasini T, Fusiello A, Trucco E, Roberto V. Making Good Features Track Better Proc . IEEE Int . Conf . on Computer Vision and Pattern Recognition (CVPR ' 98) 1998; 178- 183 and in Ruiter, N. V. ; Muller, T . O . ; Stotzka, R. ; Kaiser, W. A. Elastic registration of x- ray mammograms and three-dimensional magnetic resonance imaging data . J Digit Imaging 2001, 14 , 52-55.が開示されている。
【０００６】
異なる時間に同一対象物の2枚の平面画像を登録し、対象物の時間に対する移動量を計算することを、2次元登録と呼ぶ。同様の処理を3次元の画像セット、例えばMRIやコンピューター断層撮影（CT）スキャンの横断面画像などに用いるアルゴリズムを、３次元と呼ぶ。ランドマークや対象物の移動は、これら画像セットの中でいろんな方向に生じる。
【０００７】
比較すべき画像の異なる様式の数に依存して、アルゴリズムは単一モード登録と、マルチモード登録アルゴリズムとに分かれる。
【０００８】
２枚の連続したMRI画像を比較することは、単一モード登録である。
【０００９】
通常、単一モード登録アルゴリズムは、ほとんどの場合2次元で働く両画像内の一定のランドマーク又は複雑な形状を識別しようとする。第二ステップにおいて、両画像内の対応するランドマークの位置が比較され、それらの移動ベクトルが計算される。この情報は、第２画像のすべての画素を、新しい位置に移動させ、移動アーチファクトを除去するために用いられる。
【００１０】
「厳格な登録」は、単一ユニットとしての2次元又は３次元画像を特定の方向に移動する。「非厳格な登録」は、一定領域の単一又は複数の画素を異なる方向に移動する。
【００１１】
例えば非常にしなやかな胸の組織に関し、非厳格な登録が移動アーチファクトを除去するために必要である。
【００１２】
骨のように形成された構造の無い特別な一貫性のために、直接隣接した胸の部分は、異なる方向に移動可能である。従って、いかなる登録アルゴリズムも他の実在するアプローチも、胸の組織上のランドマークを識別し、広げることは明らかである。
【００１３】
一つのアルゴリズムグループは、有限数の発見すべきランドマークを定義する。これらは、それらの価値に従って記憶される。胸のMR画像に関し、MRIの脂肪と空気間の高いコントラストのおかげで、これらのアプローチは胸の外境界又は胸部構造部で多くの数のランドマークを見つけることがある（Huwer, S . ; Rahmel, J. ; Wangenheim, A. v . Data-Driven Registration for Lacal Deformations . Pattern Recognition Letters 1996, 11, 951-957）。
【００１４】
圧縮パネルを用いるか、又は少なくともうつ伏せの位置で横になることで両胸を固定することで、通常は圧縮装置に直接接触する肌のすべての部分を良好に固定することができる（図１）。対照的に、両胸の内部領域は、心臓の鼓動と呼吸によって僅かに動き続ける。脂肪と柔組織間の外側皮膚と空気間と比較して非常に低いコントラストにより、これらのアルゴリズムにおいて胸の中心にある有効なランドマークの数は、低く留まっている。
【００１５】
ランドマークを選択することにおける一般的な問題点は、非常に少ないランドマークの選択により、不十分な又は不正確な登録が行われることである。しかし、追加のランドマークを選択することは、正確な登録を保証するために必要ではない。なぜなら、人間の観察者又は何らかのプログラムが、例えば胸の中央の低いコントラストや複雑な組織タイプによって生じる再割り当てされた低い確実性のランドマークを含めるように強いられるからである。ランドマーク数の上昇は、常に登録計算の複雑さを著しく増加させる。
【００１６】
所謂非厳格な登録に特に適した特徴選択と特徴追跡アルゴリズムは、所謂ルーカス＆カナダアルゴリズムと、その特有のピラミッド型実行であり、詳細には以下に開示されている（Lucas BD, Kanade T . An Iterative Image Registration Technique with an Application to Stereo Vision , IJCAI81; 674-679 1981 and in Bouguet J- Y . Pyramidal Implementation of the Lucas Kanade Feature Tracker, TechReport as Part of Open Computer Vision Library Documentation, Intel Corporation）。これらは付録１として下に同封する。ルーカス＆カナダ技術と、そのピラミッド型特徴追跡実行は、特に画像中の確かなランドマークの自動識別や、異なる時間での二つの画像間のランドマークの追跡に適している。選択された特徴の置換や移動は、移動量又は所謂光学的フローベクトルとして測定される。この技術は、例えばロボット工学の画像認識などに広く応用されており、デジタル画像処理の当業者には一般的な知識の一部であり、画像認識はいくつかの基本コースで教えられている。
【００１７】
しかし、デジタル診断画像、より具体的には胸画像へのルーカス＆カナダ技術の応用の実用的な用途を考えると、この技術は、興味の外側まで含めた非常に多数の特徴を識別し、例えば胸の中の空気といったノイズの影響を受ける。
【００１８】
以下の例では、上記問題のより詳細な観点を与える。
【発明の開示】
【発明が解決しようとする課題】
【００１９】
コントラストを強調した（CE）磁気共鳴乳房撮影法（MRM）は、侵略する腫瘍に対し素晴らしい感度を有する有用で重要な調査法である。それは、胸の柔組織の密度と独立な高い負の予測値を示す。最近のマルチセンター研究では、感度と特異性は、画像解析技術に従って異なることが示されている。1.5Tスキャナーを用いた、９６、９３又は８６％の感度で３０、５０又は７０％の特異性のものがそれぞれ計算されている。従って。MRMの主な欠点は、誤った好結果の診断が高い確率で生じることである。
【００２０】
通常MRイメージングは、1.0又は1.5テスラ撮像装置でメーカーの二重ブレストコイルを用いて（米国では、しばしば単一ブレストコイルが用いられる）、患者をうつ伏せにして行われる。プロトコルは、少なくともあらゆる向きで捕捉時間３０〜１８０秒でのダイナミックコントラストで強調されT1で重み付けされたシーケンスを有する。少なくとも一つの測定が、ガドペンテト酸ジメグルミンの大丸薬注入又は他の常磁性のMRコントラスト薬剤の投入前に行われ、幾つかの測定がその後に行われる。
【００２１】
胸の組織の領域への実際のコントラスト薬剤の取り込みを計算するために、コントラスト薬剤投与前の各画素の信号強度は、コントラスト薬剤投与後の信号強度から引き算されなければならない。MR画像のスライス厚さで生じる部分的な体積効果と同様に、呼吸又は心臓の鼓動によって生じる両胸の最小限の移動によって、異なる時間に得られた同じ画像位置の画素は、正確に胸の組織の同じ厚さを示さない。従って引き算された画像は、絶対的に黒にはならない（腫瘍部を除く）。少しの移動の効果は、胸の外側境界を良く説明することができるであろう。脂肪組織の高い信号強度と周囲の空気のおよそゼロの信号強度とによって、非常に小さな移動が、脂肪を空気の以前の画像位置に置く。結果として、引き算の画像は、画像化された胸の少なくとも部分周りに取り込まれた強いコントラスト薬剤を現わす厚い白い線を示す（図２で、胸の白い外側境界がサインによってマークされている）。
【００２２】
デカルト体系によって定義された3次元又は体積画像を得ることを考慮すると、ｘ、ｙ方向へのあらゆる移動に加えて、（胸に向かう）ｚ方向へのいくらかの移動が常に見られる。これは、一連の調査中の患者のくつろぎによって生じる。
【００２３】
アーチファクトを減らすために、胸は通常ある種の圧縮法（メーカーに依存する）によって胸コイルの内側に固定される。しかし、この事実にかかわらず、最小の移動は常に発生する。従って、胸を固定しているにもかかわらず、画像には僅かな移動アーチファクトがｘ、ｙ、ｚ方向に見られる。移動アーチファクトは明るい画素として引き算画像の中で視覚化される。移動アーチファクトがなければ、引き算画像は腫瘍組織に占領された領域を除き、絶対的に黒である。
【００２４】
一定数のランドマークが胸の外側の周囲の空気のノイズの大きな信号の中から検出される。考慮された各特徴が追跡されるべきであり、不必要な計算上の負荷が発生するという事実の他に、ノイズはランダムな効果であって、ノイズに含まれる全ての特徴又はランドマークは、第２の、又は後に続く画像の対応するランドマークを有さない。しかし事実として、ノイズ内で見出されたランドマークに対応するアルゴリズムがサーチする第２画像やそれに続く画像内のランドマークは、ともかく第２画像や画像セット内に一定量のランドマークを再割り当てする。これが、画像の間違った登録に繋がり、結果に影響を与える。
【課題を解決するための手段】
【００２５】
本発明の目的は、画像、特に3次元生物医学的診断用画像を登録する方法に存する。この方法は、3次元での移動追跡を可能とし、デジタル画像内の特に興味のある領域の外側のノイズによるアーチファクトを除去し、隣接する組織の信号強度の高低による差異と独立に胸にかかるすべてのランドマークを広げることによって、既知の登録方法の欠点を効果的に克服するものである。
【００２６】
本発明によれば、画像を登録する方法は以下のステップからなる。
ａ）画素の2次元又は３次元アレイによって形成された画像からなる、解剖学上の部分の組織又は組織帯をＭＲＩ、超音波、又は放射線で走査して得られた第１デジタル画像又は横断面画像セットを提供するステップと、
第二の時間に同じ解剖学上の部分の組織又は組織帯をＭＲＩ、超音波、又は放射線で走査して得られた少なくとも第二デジタル画像又は横断面画像セットを提供するステップと、
ｂ）前記第１画像又は画像セットの中で、一定数のランドマーク又は所謂特徴を、一定数の画素を選択することによって定義し、追跡すべき特徴のリストを生成するステップと、
ｃ）前記第１又は第２画像又は画像セットの中から特徴として選択された画素を有する画像の光学的フローベクターが定義されることによって、特徴として選択された各画素の位置が追跡されるステップと、
ｄ）逆光学フローを第二画像又は画像セットに加えることによって、第一と第二画像又は画像セットが登録されるステップと、
ｅ）ステップｄ）での登録後に、第２画像の第２画像データから第１画像の画像データを引き算するステップと、
ｆ）第２画像の第２画像データから第１画像の画像データを引き算して得られた画像データを表示するステップ。
【００２７】
ステップｂ）の特徴選択は、以下により詳細に記載する。
Ｂ１）第一画像又は第一断面画像セットの各画素周りの隣接する画素領域が定義され、各画素領域が制限された数の画素を有し、単一画像の場合、２次元領域が選択され、断面画像セットの場合、３次元領域が選択されるステップと、
Ｂ２）各画素に対して、前記画素領域の全画素の平均信号強度が求められるステップと、
Ｂ３）平均信号強度閾値が求められるステップと、
Ｂ４）ステップ２で求められた各画素領域の平均信号強度と、前記平均信号強度閾値とが比較されるステップと、
Ｂ５）前記領域の平均信号強度が、ステップ４の閾値より大きい場合、画素は追跡されるべきランドマークとして定義され、追跡されるべきランドマークリストに追加されるステップとから成る自動特徴選択ステップ。
【００２８】
平均信号強度閾値は、調整可能であり、変更可能である。かかる閾値を決めるために、経験的な基準を適用することができる。例えばテスト画像又はサンプル画像に本方法を実行し、異なる閾値に対して得られた結果を比較する。
【００２９】
本発明によれば、画像の追跡すべき特徴を定義できるすべての種類の画素を選択する方法を適用可能であるが、最良の結果は既知のルーカス＆カナダ法による自動特徴選択法によって得られ、この方法の詳細は、Lucas BD, Kanade T . An Iterative Image Registration Technique with an Application to Stereo Vision, IJCAI81; 674-679 1981に記載されている。
【００３０】
所謂光学的フローという移動又は置換ベクトルによって、第１画像から第２画像への特徴の移動を追跡することは、既知の方法である。またこの場合、本発明をランドマークや特徴を利用する幾つかの追跡方法に適用することができる。特に、もし診断画像撮影中に制限された胸の小さな移動に関わるものである場合、そしてルーカス＆カナダ特徴選択法を使用する場合、本発明が、好ましくは所謂「ルーカス＆カナダ特徴追跡のピラミッド型実行」（付録１参照）を使用する。この「実行」は、当業者には既知であり、幾つかの大学の教育課程で教えられている。このルーカス＆カナダ特徴追跡のピラミッド型実行の特定の、詳細な内容は、Bouguet J-Y . Pyramidal Implementation of the Lucas Kanade Feature Tracker, Tech Report as Part of Open Computer Vision Library Documentation, Intel Corporation and in Shi J, Tomasi C . Good Features to Track . 1994 IEEE Conference on Computer Vision and Pattern Recognition (CVPR' 94) 1994; 593- 600 and in Tommasini T, Fusiello A, Trucco E, Roberto V. Making Good Features Track Better . Proc . IEEE Int . Conf . on Computer Vision and Pattern Recognition (CVPR ' 98) 1998; 178-183.に開示されている。
【００３１】
ひとたび全ての特徴に対して光学的フローベクトルが定義されると、逆フローを第２画像に加えることで、第２画像の各特徴とその周辺画素とを第１画像の特徴位置に移動することによって、２つの画像の登録が実行される。
【００３２】
既に述べたように、上記方法は、磁気共鳴イメージング（MRI）又はコンピュータ断層撮影法（CT）の同じ又は異なる種類のスキャナーを用いた、２次元又は３次元画像を登録するために用いることができる。
【００３３】
上述の登録方法は、MRI又は超音波画像といったコントラスト媒体強調診断イメージング、特にコントラスト強調磁気共鳴乳房X線撮影の方法を組み合わせることでも提供される。この方法において、胸の組織の領域への実際のコントラスト薬剤取り込みを計算するために、コントラスト薬剤塗布前の各画素の信号強度が、コントラスト薬剤塗布後の信号強度から引き算されなければならない。このように、画像化された組織においてコントラスト薬剤が存在しない時間における第１画像と、コントラスト薬剤が存在する第２の又は後の時間の第２画像とを登録することと、画像化された組織内を一面に覆うことが、もし登録が実行されなかった場合、患者の移動によって画像データの引き算に生じるアーチファクトを、高い度合いで抑圧するために役立つ。主となる考えは、２枚の画像間を引き算することで、コントラスト薬剤の溜まった組織に対応した画素が高い信号強度を有することである。コントラスト薬剤の存在しないその他の部分は、低い信号強度の画素として、特に暗い灰色又は黒色で表示される。
【発明を実施するための最良の形態】
【００３４】
本発明の一つの観点は、ルーカス＆カナダアルゴリズムの２次元から３次元への拡張であるが、以下のグラフィック例は、理解を容易にするために２次元画像を用いている。３次元画像に対する本方法の一般化は、数学的な形式論を用いて概要によって与えられる。
【００３５】
３次元と３次元データを現す断面画像セットの自動特徴選択作業は、ルーカス＆カナダの２次元特徴検出アルゴリズム実体的な拡張であって、以下に詳細に説明する。
【００３６】
時間を少し変えて記録した同じ患者の２枚の画像を、IとJと定義する。
【００３７】
最初の画像Iにおいて、２枚の画像間で追跡すべき特徴であると考えられる一定数の選択された画素が識別されなければならない。
【００３８】
Iをオリジナル／最初の画像とし、Jを次の画像とし、ここで対応する特徴は見出されていなければならない。
【００３９】
第一のステップは、最初の画像の最小固有値の計算からなる。
【００４０】
これは、以下のように実行される。
第１画像を、３つの勾配量Ix、Iy、Izに変換し、（デカルト座標系のｘ、ｙ、ｚ方向で画像を描写する）、一方ｘ、ｙ、ｚは各単一画素の座標である。
勾配量は、隣接する画素に対する各画素の強度勾配から算出され、隣接する画素は一般的にx方向で（２？ｘ＋１）のサイズ、ｙ方向で（２？ｙ＋１）のサイズ、ｚ方向で（２？ｚ＋１）のサイズを有する所謂ウインドウからなる。ここで、？ｘ、？ｙ、？ｚ＝１、２、３、４、．．．、ｎ 画素である。
【００４１】
目標画素の中央の３×３×３アレイの隣接画素のサイズを考慮して、３つの勾配量は以下のように定義される。
ａ）ｘ方向の勾配量
ｂ）ｙ方向の勾配量
ｃ）ｚ方向の勾配量
各画素の次のステップとして、傾斜マトリックスＧが、上述の３つの計算された勾配量を用いて表される。
傾斜マトリックスは、以下のように生成される。
各画素の各勾配マトリックスに対し、最小固有値？minが、以下に従って傾斜マトリックスより算出される。
勾配マトリックスGの固有値は、以下のように計算される。
閾値は、固有値？ｍの最小値として定義される。
画素が追跡可能な特徴を現しているか否かを考慮するために、以下の基準が適用される。
１．もし？ｍが閾値より大きい場合、画素の実際の位置が追跡すべき特徴リストに保存される。
２．もし？ｍが固有値？ｍの全ての最小値に対する最大値のパーセンテージより小さい場合、追跡すべき特徴リストから落とされる。
３．実際の特徴の周りの定義された領域（調整可能な、例えば３×３×３）に他の特徴が存在する場合、固有値？ｍのより小さい最小値を有する特徴は追跡すべき特徴リストから落とされる。
４．もし特徴位置周りの３ｄ隣接領域（例えば？＝１）の平均信号強度が、ゼロに近い非常に低い値の（例えば１０）平均信号強度閾値より小さいか又は等しい場合、この特徴は、胸の外側の周辺の空気に位置すると見なされ、追跡すべき特徴リストから落とされる。
５．最後に、定義可能な最大数の特徴が保持される。もし追跡すべき特徴リストに、それ以上の特徴があった場合、より小さい最小値？ｍを有する特徴がリストから落とされる。
【００４２】
これらのステップは、２次元の場合を参照して図３ないし６で視覚的に説明する。この場合、画素アレイは小さな四角のアレイとして示されている。異なる画素は、第一基準を満足する固有値の最小値を有する選択された特徴の位置に対応して示される。これらの特徴と位置は、対応する画素に異なる様相を与えることによって、図３で強調表示されている。図４において、５×５の隣接領域が、その内側に円を書くことで強調表示されている。
【００４３】
図１において、患者の胸を固定する典型的な様相が図示されている。通常、MRイメージングは、患者がうつ伏せ状態で、メーカーの２重胸コイルを用いて（しばしば米国においては単一の胸コイル装置が用いられる）、1.0から1.5テスライメージャーを用いて行われる。プロトコルは、少なくともダイナミックコントラスト強調T1重み付けシーケンスをすべての方向に対して利得時間３０〜１８０秒で有する（個々のプロトコルで異なり、脂肪抑圧を付加的に使用可能である）。少なくとも一つの測定が、0.1mmol/Kgから0.2 mmol/kg体重のガドペンテト酸ジメグルミンの大丸薬投与又は他の常磁性MRコントラスト薬剤投与前に、幾つかの測定が投与後に行われる。
【００４４】
少なくとも第一MRI画像は、コントラスト媒体が肘前の静脈に投与され、胸組織に達する前の、時間t0の時に得られる。少なくとももう一つの画像が、コントラスト媒体投与後の時間t1に得られる。通常、時間を遅らせた一連の画像が得られる（例えば６回以上連続して同じ対象に対しコントラスト媒体投与後に得る）。
【００４５】
現在の画像を単純化するために、断面図の3次元画像から取り出した２枚の画像は、一定の(スライス)厚さを有する。従って、2次元画素に関して論じる代わりに、3次元画素に関してより詳細に論じる。各画素は、対応する信号強度と相関関係のある所定のグレーレベルを有する。図３ないし１０の例において、グリッドは画素アレイを表現している。グリッドの各要素は、画像を形成する画素である。
【００４６】
図３において、例として４種類の画素がランドマークとして特徴づけられており、対応する勾配マトリックスの最小固有値に従って保存される。102と示された画素は、最大固有値を有する画素であり、202と示された画素は、大きな固有値を有する画素であり、302は小さい固有値を有する画素を、402は最小固有値を有する画素を示す。
【００４７】
図４において、各画素に対して直径５画素のEucliadian距離内の隣接領域が特徴として選択された中央画素の周りに定義される。この隣接画素の選択は、円として図４に示されている。
【００４８】
図４において、このイベントの第一の場合が502'で示された厚い円で強調表示されている。ここで、特徴として選択された２画素は、画素のEucliadian距離として上で定義されたと同じ隣接領域内にある。これらの2画素は、202'と302'で示され、上記定義に従って、特徴202'として選択された画素は、第二の画素302'より大きな固有値の勾配マトリックスを有する。従って、この第二画素302'は除かれる。
【００４９】
図５は次のステップを示し、図３で特徴として選択された302'はもはや存在しない。図５において、画素202'はその隣接領域を示す円を備えている。図４と同様に、202''の隣接領域に202''より小さな固有値を有する傾斜マトリックスの画素302''が存在する。
【００５０】
従って、この場合も特徴として選択された画素302''は追跡すべき特徴のリストから除かれる。
【００５１】
図６の画素アレイにおいて、少なくとも幾つかの特徴として選択された当初からの画素が存在する。単一の隣接領域の中では、一つ以上の特徴はもはや見出すことができない。従って追跡すべき特徴は、画素全体に広がった。
【００５２】
本発明の更なるステップによれば、画素の強度は更に調査され、所定値より低い強度を有する特徴として選択された画素もまた、追跡すべき特徴のリストから除かれる。ほぼゼロの非常に低い閾値が、周辺の空気内で見つけられる特徴の特徴づけに用いられる。このステップは図示されないが、当業者には自明である。
【００５３】
ひとたび画素が2次元または３次元画像の特徴として考慮すべきと決定すると、特徴の追跡が実行される。好ましき実施例によれば、これは所謂「ルーカス・カナデ特徴追跡のピラミッド型実行」を用いて実行される。詳細は、Jean Yves Bouguet in the article "Pyramidal Implementation of the Lucas Kanade Feature Tracker, Description of the Algorithm" published By Intel Corporation, Microprocessor Research Labsに記載されている。
【００５４】
上述の出版物において2次元例として開示されている理論を、以降3次元の場合に更に適用する。以下のステップは、上述の方法とアルゴリズムに従って個性化された追跡すべき特徴を現す各画素に対して実行されなければならない。
【００５５】
ボリュームIのポイントｕと、ボリュームＪの対応するポイントｖを定義する。第一ステップにおいて、二つのピラミッドが画像ＩとＪに対しレベルＬ＝０、．．．、ｍを現すＩLとＪLとして、また当初ボリュームを表現する基部（Ｌ＝０）を有して組み立てられなければならない。各階のボリュームは、直接隣接する領域の一定量の画素を一つの平均値と組み合わせることによってサイズが減少する。追加の階の数は、ボリュームのサイズに依存する。続くステップにおいて、所謂グローバルピラミッド型移動ベクトルｇが最も高い階で初期化される。
【００５６】
全てのレベルＬに対して、結果としての移動ベクトルｄLが計算される。
【００５７】
これは以下のステップに従って、実行される。
ａ）各ボリュームｌLにおいて、ポイントｕLの位置が以下の式で計算される。
ここでｐは実際のボリューム位置である。
【００５８】
ｂ）ボリュームのピラミッド型の表現の各レベルＬで、勾配はデカルト座標の各方向ｘ、ｙ、ｚに対して計算される。
【００５９】
デカルト座標ｘに対し、ボリューム勾配が下式に従って計算される。
【００６０】
この式は、他の二つのデカルト座標ｙとｚに対しても同様に適用され、
【００６１】
ｃ）ボリューム勾配は、次に以下の式に従って勾配マトリックスを計算するために用いられる。
【００６２】
ここで、値？は、画素に影響を与える隣接領域のエリアサイズを定義する。
【００６３】
マトリックスの要素は上述の定義に従った方法と同様に、下記式で得られる。
【００６４】
ｄ）各レベルＬに対して、相互作用ベクトル？が初期化される
【００６５】
【００６６】
以下の式に従って、ミスマッチベクトルが計算される。
【００６７】
光学的フローベクトルは、下記式によって定義される。
【００６８】
相互作用移動ベクトルは上で定義した光学的フローベクトルを用いて下のように計算できる。
【００６９】
そして、これがレベルＬに対する最終光学的フローと等しいセットである。
【００７０】
上のステップは、最終レベルＬ＝０まで各レベルで繰り返される。
【００７１】
次のレベルＬ−１に対し、グローバル光学フローは上記計算値から計算できる。
【００７２】
ｆ）最終光学的フローベクトルは、次式となる。
【００７３】
ｇ）ボリュームＩのポイントＵに対応するボリュームＪのポイントｖの座標は、下記のように計算できる。
ｖ＝ｕ＋ｄ
【００７４】
図７ないし１０は、現実の表現を明らかに単純化した上述の追跡方法のグラフィカルで模式的なアイディアを与えるものである。
【００７５】
図７ないし１０において、異なる時間Ｔ０とＴ１に撮影された二つの画像ＩとＪが、図３ないし６の例と同様に、画素の2次元アレイとして現されている。
【００７６】
上述の方法に従って追跡されるべく選択された特徴に対応する画素２０が、画像Ｉの２次元画素アレイ内に示されている。同じ位置の画像Ｊにおいて、画素２０’’は画像化された対象、この場合胸が動かない場合に画素２０と対応する。画素２０’’の隣接する特定の画素が、中央が画素２０’の７×７の画素アレイに対応してグレイの影付き四角で示されており、境界部は番号１０で示されている。
【００７７】
２０’は、画像Ｉの特徴２０の新たな位置を示している。これは、画素２０の位置にある胸の組織が時間と共に画素２０’の位置に動いたことを意味する。明らかに、画素２０’は追跡が実行されたことを示すためにのみ図示されている。実際の場合、この画素は知られず、又は追跡が不要である。
【００７８】
数学的な形式論によって既に述べられた選択された特徴に対応する画素を追跡する一般的な方法は、以下に簡単のために2次元の場合に限定して図面とともに実際的な例を記載する。しかし、当業者は追跡の一般的な数学的記載がそれでも理解でき、3次元の場合に追跡する方法を評価する。
【００７９】
追跡の初めに、画像Ｉの画素２０は、同じ位置、即ち画像Ｊの画素２０’’にあるものと仮定する。この画素２０’’は以降、被追跡画素と呼び、特徴として選択された一方画像Ｉの画素２０は特徴画素と呼ぶこととする。第一に、被追跡画素２０’’の隣接領域に対応する一定領域の逆勾配マトリックスにアルゴリズムを適用する。本例では画像Ｊの画素２０’’を中心とした７×７画素アレイ（この値は調整可能）を計算する。所謂追跡領域にこの隣接領域は対応する。ルーカス＆カナデの特徴追跡アルゴリズムのピラミッド型実行を適用することによって、3次元画像にまで適用でき、その詳細は後述するが、画像Ｊの追跡された特徴の新しい位置を見つけることができる。
【００８０】
画像Ｊにおける画素２０’’の追跡される特徴の新しい位置は、光学的フローベクトルを決定する。これは、第一の反復ステップで計算されることによって定義される。図７と図８の比較から明らかなように、追跡領域は、同じ組織領域２０’の新しい位置により近い追跡領域の中央になるように移動する。
【００８１】
追跡は、画像Ｉの特徴を現す画素２０と、画像Ｊの追跡領域中央の新たに識別された位置との間でミスマッチベクトルの計算を繰り返すことによって、さらに反復的に適用される。これが新しい光学的フローベクトルの定義をもたらし、図９に示すように画素隣接領域又は追跡領域の新たな移動をもたらす。
【００８２】
画像Ｉの特徴画素２０と選択領域の中央との間のミスマッチベクトルを決定するステップは、ミスマッチベクトルが調整可能な最小値に達するまで行われる。
【００８３】
この状態は、図１０に示される。図７ないし９の図で表現されているものに関し、画素２０’はその位置を維持しており、追跡領域が移動しているという事実を強調している。追跡領域は、最終的に正しい新しい位置２０’を中心とした位置に移動し、画像Ｉの画素２０としての胸の同じ位置を表現している。
【００８４】
上述の例は、ただ一つの特徴が画素によって表現されている場合に実行される。上述のステップは、画像内で選択された各特徴に対して実行されなければならない。
【００８５】
全ての特徴がその新しい位置まで追跡された後、画像上の全ての異なる位置にある、胸の組織を表現する領域内にある一定量の光学的フローベクトルが、存在する。他の登録アルゴリズムと比較して最大の長所は、このアルゴリズムが画像全体に光学的フローベクトルを広げることである。
【００８６】
次に、逆光学的フローベクトルを画像Ｊの各特徴とその周囲の隣接領域とに加えることで、モーフィングを実行する。これによって、新しい視覚的画像である画像Ｊ’となる。
【００８７】
図１１は、この方法の実際の診断画像での効果を説明している。左側の第一の画像セットは、コントラスト薬剤投与無しでの時間t0での胸の代表的なスライスを示す。左側の第二画像セットは、コントラスト薬剤投与後の時間t1での同じスライス部位を示す。左側の第三画像セットは、時間t1での各画素の信号強度から、時間t0での強度を引き算した、所謂引き算画像の結果を視覚化したものである。時間間隔t1−t0の間に相対的な胸の移動の無い場合は、結果画像は画素の広い部分でコントラスト薬剤強調の生じない黒で示される信号強度ゼロの部分が含まれる。時間間隔t1−t0間の移動は、組織の画像の同じ領域とは対応しない画像ＩとＪの同じ位置にアーチファクトが生成されるという効果を得る。本例において、引き算画像は小さな矢印で示されたコントラスト強調された腫瘍を含む一つの胸を示している。ほとんどの種類の胸の悪性と良性の腫瘍は、血管分布の増加によって生じる追加的なコントラスト薬剤の取り込みを示す。両方の胸の上部の追加的な白い構造体は、移動アーチファクトによるものである。両方の胸の中央に残った白い構造体は、移動との組み合わせによる通常の胸の柔組織の弱いコントラスト薬剤の取り込みの結果である。
【００８８】
図１１の右側は、登録アルゴリズムが適用された対応する画像を示している。
【００８９】
時間t0での画像は、アルゴリズムによっては変更されないため、同じままである。しかし、3次元特徴追跡後は、新たな視覚的画像が計算される。時間t1'での右側の画像は、登録適用後の対応する同じスライス位置のバーチャル画像を示す。右側の最後の画像は、時間t1'での各画素の信号強度から時間t0での信号強度を引き算した、画素対画素の引き算結果を新たな（バーチャル）引き算画像として視覚化したものである。一つの胸は、小さな矢で示すコントラスト強調腫瘍を示している。通常の胸の柔組織を現す胸の範囲内の領域では、いくつかの中間のコントラスト薬剤取り込みが見られ、移動アーチファクトの減少により少し暗く見えている。
【００９０】
図１１の例は、コントラスト薬剤投与前後の2枚の連続画像のみにより簡略化したものである。
【００９１】
図２から図１１の画像はすべて、同じMRスキャナで撮影し、登録以外の加工は行っていない。
【００９２】
実際には、連続した断面写真（例えば50枚以上）がコントラスト媒体の注入前後の通常の時間間隔で得られ、引き算画像がそれぞれの連続画像の組み合わせに対して計算される。
【００９３】
一連の画像の中で一枚以上の横断面画像を得るのに加えて、2連以上（例えば５以上）の画像の中から次々に得ることができる。これは、幾つかの引き算画像となる。もちろん、開示したアルゴリズムは、コントラスト薬剤移動の前後に関わらず、連続写真の全ての組み合わせに適応でき、その間に生じたいかなる移動も減じることができる。
【００９４】
以下の別表において、ルーカス＆カナデの追跡アルゴリズムのピラミッド型実行の詳細な理論を示す。
【００９５】
本発明に係る上述の方法は、コンピューターで実行可能なプログラムとしてコード化されたアルゴリズムとして実行される。
【００９６】
このプログラムは、例えば磁気的、光磁気的、又は光学的基板、例えば1枚以上のフレキシブルディスク、1枚以上の光磁気ディスク、1枚以上のCD-ROM又はDVD-ROM、その他代替的なものなどの携帯型又は静止型記憶装置や、携帯型又は据え置き型ハードディスクなどに保存可能である。
【００９７】
擬似コードで書かれたコード化されたプログラム例を以下に示す。
【図面の簡単な説明】
【００９８】
【図１】胸部MRIで用いられる胸固定手段の概要図である。
【図２】引き算過程と、登録ステップが実行されない場合の移動アーチファクトの発生を示す図である。左上と右上の図は、同じMRスキャナ位置で撮った２枚の胸の断面図である。第２画像は、第１画像撮影の後、約２．５分のコントラスト媒体投与後に測定した。両画像は、断面画像の３次元セットの代表部分である。
【図３】模式的な例と共に特徴選択基準が示された模式的な２次元画像である。
【図４】模式的な例と共に特徴選択基準が示された模式的な２次元画像である。
【図５】模式的な例と共に特徴選択基準が示された模式的な２次元画像である。
【図６】模式的な例と共に特徴選択基準が示された模式的な２次元画像である。
【図７】どのように特徴追跡が実行されるかを模式的に示した図である。
【図８】どのように特徴追跡が実行されるかを模式的に示した図である。
【図９】どのように特徴追跡が実行されるかを模式的に示した図である。
【図１０】どのように特徴追跡が実行されるかを模式的に示した図である。
【図１１】コントラスト媒体投与前後の異なる時間に撮影した胸の２枚の画像と、特徴追跡と登録実行後の引き算をした画像であり、左側には比較として図２の画像例を示す。
【特許請求の範囲】
【請求項１】
ａ）画素の2次元又は３次元アレイによって形成された画像からなる、同一対象物の少なくとも第１と第２デジタル画像又は横断面画像セットを提供するステップと、
ｂ）前記第１画像又は画像セットの中で、一定数のランドマーク又は所謂特徴を、一定数の画素を選択することによって定義し、追跡すべき特徴のリストを生成するステップと、
ｃ）前記第１又は第２画像又は画像セットの中から特徴として選択された画素を有する画像の光学的フローベクターが定義されることによって、特徴として選択された各画素の位置が追跡されるステップと、
ｄ）逆光学フローを第二画像又は画像セットに加えることによって、第一と第二画像又は画像セットが登録されるステップとからなる、
対象物の移動によって画像アーチファクトを減少させるための生物医学的画像の登録方法であって、
Ｂ１）第一画像又は第一断面画像セットの各画素周りの隣接する画素領域が定義され、各画素領域が制限された数の画素を有し、単一画像の場合、２次元領域が選択され、断面画像セットの場合、３次元領域が選択されるステップと、
Ｂ２）各画素に対して、前記画素領域の全画素の平均信号強度が求められるステップと、
Ｂ３）平均信号強度閾値が求められるステップと、
Ｂ４）ステップ２で求められた各画素領域の平均信号強度と、前記平均信号強度閾値とが比較されるステップと、
Ｂ５）前記領域の平均信号強度が、ステップ４の閾値より大きい場合、画素は追跡されるべきランドマークとして定義され、追跡されるべきランドマークリストに追加されるステップとから成る自動特徴選択ステップが実行されることを特徴とする方法。
【請求項２】
前記平均信号強度閾値が実験によって決定され、試験セット又はサンプル画像に対し前記方法を実行し、得られた異なる閾値の結果を比較することによって決定されることを特徴とする請求項１記載の方法。
【請求項３】
第２画像内の追跡されるべきランドマークに選択された画像の画素が、所謂ルーカス＆カナダ自動特徴追跡方法又はアルゴリズムによって自動的に追跡されることを特徴とする請求項１又は２記載の方法。
【請求項４】
請求項３記載の２次元画像にのみ導入されたルーカス＆カナダアルゴリズムの初期の実行が、２次元画像に代わって３次元画像にまで拡大され、
ａ）３次元画素により形成された第１の３次元画像ボリュームＩが提供されるステップと
、
ｂ）各目標画素に対して各画素を取り囲み、調整可能な例えば前記目標画素を中心に３×３×３アレイの画素を有する画素領域が定義されるステップと、
ｃ）デカルト座標系を定義し、前記デカルト座標系の各画素に対し相対的に、隣接画素に対して相対的に第１画像ボリュームの各画素の強度傾斜によって形成される下記式で定義される所謂傾斜ボリュームが計算されるステップと、
ｄ）各目標画素に対する以下の式で定義される所謂傾斜マトリックスが計算されるステップと、
ｅ）３次元画像ボリュームの各目標画素の各傾斜マトリックスに対し、以下の式を適用することで、最小固有値が計算されるステップと、
定義
傾斜マトリックスの最小固有値の計算
ｆ）最小固有値の閾値が定義されるステップと、
ｇ）対応する傾斜マトリックスが、
ｉ）固有値が閾値より大きい
ii）固有値が固有値のすべての最小値に対する最大値の割合より小さくないこと
iii）もし特徴として他の画素が、定義された選択された特徴としてまた維持された目標画素周りの画素領域内に存在する場合、傾斜マトリックスがより大きな最小固有値を有する画素のうちの一つのみが、特徴として選択され、他の一つは、追跡可能な特徴のリストから破棄される
という条件を満足する最小固有値を有する画像Ｉの追跡すべき特徴としての各画素が選択されるステップと、からなる生物医学的横断面画像ボリュームに適用されることを特徴とする方法。
【請求項５】
特徴として選択された画素の隣接領域の平均信号強度に対する閾値が約１０にセットされていることを特徴とする請求項４記載の方法。
【請求項６】
第１画像ボリュームＩの特徴として選択された各画素に対し、前記特徴が第２の遅い時間に同じ対象物の別の画像Ｊの画素アレイ内の位置に対して相対的に追跡され、前記追跡は、所謂ルーカス＆カナダの特徴追跡のピラミッド型実行によって実施されることを特徴とする請求項１ないし５のいずれか1記載の方法。
【請求項７】
画像Ｉの３次元画素アレイに対応する画像ボリュームＩ内のポイントとして、ポイント
ｕを、画像Ｊの3次元画素アレイに対応する画像ボリュームＪの対応するポイントをｖと定義し、これらのポイントが画像Ｉの特徴として選択された画素の座標を有するステップ
と、
当初のボリュームを表す基部（Ｌ＝０）と共に、レベルＬ＝０、．．．ｍまでのボリュームを表す、画像ボリュームＩとＪの２つの理想的なピラミッドＩLとＪLとが構築されるステップと、
直接隣接する領域の画素の一定量を一つの平均値と組み合わせることで、各追従する階のボリュームのサイズが減少されるステップと、
もっとも高い階Ｌmで初期化される所謂包括的ピラミッド型移動ベクトルｇを定義するステップと、
ｉ）ｕLとしてのポイントｕの位置を定義し、以下の式で前記位置を計算するステップ
と、ここでＬはピラミッドの実際の階又はレベルでの値であり、ここでｐは実際のボリューム位置であり、
ii）ボリューム傾斜のための以下の式に従って、デカルト座標ｘ、ｙ、ｚの各方向に対
してボリューム傾斜を計算するステップと、
iii）以下の式に従って傾斜マトリックスを計算するために傾斜ボリュームを用いるステップと、
ここでの値は、追跡された特徴を表現する目標画素に影響を与えるエリアサイズ又は画素の隣接領域を定義し、
iv）各レベルＬに対し、下の式で定義される反復ベクトルを初期化するステップと、
ｖ）kが１から最大カウント数Ｋまで、又は最小置換
まで、以下のボリューム差を計算するステップと、
vi）以下の式に従って、ミスマッチベクトルを計算するステップと、
vii）以下の式に従って、光学フローベクトルを定義するステップと、
ここでＧ-1はステップiii）で定義された逆傾斜マトリックスＧである。
viii）上で定義された光学フローベクトルを用いた双方向ベクトルを計算する以下のステップと、
これをレベルＬの最終光学フローと等しくする以下のステップと、
ix）前記ステップｉ）からviii）を各レベルで最終レベルＬ＝０になるまで繰り返し、下記式で定義される最終光学フローベクトルに達するステップと、
ｘ）追跡されるべき特徴を現すポイントｕに対応するボリュームＪのポイントｖの座標
を下の式によって定義するステップと、
xi）画像Ｉ内で選択された特徴に対応する各画素に対して、上記方法を繰り返すステ
ップと、
からなる所謂ルーカス＆カナダの特徴追跡のピラミッド型実行であることを特徴とする請求項６記載の方法。
【請求項８】
逆光学フローベクトルを画像Ｉの選択された特徴に対応する画素のポイントｕに対応する特徴を現す画素に対応して識別される画像Ｊの各ポイントｖに加えることによって、 又は光学フローベクトルをそれぞれ画像Ｉの選択された特徴として識別される画素に対応したポイントｕに加えることによって、２つの画像ＩとＪの登録が実行されることを特徴とする請求項７記載の方法。
【請求項９】
造影剤強化診断イメージングのための方法との組み合わせ、特に造影剤強化ＭＲＩ又は超音波イメージングを提供することを特徴とする方法であって、
ａ）ＭＲＩ又は超音波撮影で得られた、少なくとも第１と第２のデジタル画像又は同一対象物の断面画像セットを提供するステップであって、
前記画像が画素の2次元又は３次元アレイによって形成され、
同一組織又は同一組織帯又は同一の解剖学上の部分の走査が造影剤の存在下で第２の時間又はいずれかの以後の時間に実行され、
ｂ）第１画像又は画像セットの内部で、特定数のランドマークや所謂特徴を特定数の画素を選択することによって定義し、追跡すべき前記特徴のリストを生成するステップと、
ｃ）特徴として選択された各画素に対し、第１画像から第２画像まで又は画像セットから前記光学フローベクトルを定義することによって、第１画像から第２画像で特徴として選択された各画素又は画像セットの位置を追跡するステップと、
ｄ）前記逆光学フローを第２画像の画素又は画像セットに加えることで、第１又は第２画像又は画像セットを登録するステップと
からなる請求項１ないし８のいずれか１に記載の方法。
【請求項１０】
Ｂ１）第１画像又は第一断面画像セットの各画素周りの画素近傍領域を定義するステップと、
ここで前記画素近傍領域は、限定された数の画素で構成され、
２次元近傍領域が、単一画像の場合に選択され、
３次元近傍領域が、断面画像セットの場合に選択され、
Ｂ２）各画素に対し、前記画素の近傍領域の全画素の平均信号強度を定義するステップと、
Ｂ３）平均信号強度閾値を定義するステップと、
Ｂ４）ステップＢ２で求めた平均信号強度と、前記平均信号強度閾値を比較するステップと、
Ｂ５）前記近傍領域の平均信号強度がステップ４で求めた閾値より高い場合、画素が追跡されるべき特徴として定義され、追跡されるべき特徴のリストに追加されるステップと、
からなる更なる自動特徴選択ステップが、特徴選択ステップの後、特徴追跡ステップの前に実行されることを特徴とする請求項１ないし９のいずれか１に記載の方法。
【請求項１１】
２つのデジタル画像又は画像ボリュームを登録するためのコンピュータプログラムであって、請求項１ないし１０のいずれか１に記載の方法がコード化され、以下の擬似コードで特徴付けられるコンピュータプログラム。
【請求項１】
ａ）画素の2次元又は３次元アレイによって形成された画像からなる、同一対象物の少なくとも第１と第２デジタル画像又は横断面画像セットを提供するステップと、
ｂ）前記第１画像又は画像セットの中で、一定数のランドマーク又は所謂特徴を、一定数の画素を選択することによって定義し、追跡すべき特徴のリストを生成するステップと、
ｃ）前記第１又は第２画像又は画像セットの中から特徴として選択された画素を有する画像の光学的フローベクターが定義されることによって、特徴として選択された各画素の位置が追跡されるステップと、
ｄ）逆光学フローを第二画像又は画像セットに加えることによって、第一と第二画像又は画像セットが登録されるステップとからなる、
対象物の移動によって画像アーチファクトを減少させるための生物医学的画像の登録方法であって、
Ｂ１）第一画像又は第一断面画像セットの各画素周りの隣接する画素領域が定義され、各画素領域が制限された数の画素を有し、単一画像の場合、２次元領域が選択され、断面画像セットの場合、３次元領域が選択されるステップと、
Ｂ２）各画素に対して、前記画素領域の全画素の平均信号強度が求められるステップと、
Ｂ３）平均信号強度閾値が求められるステップと、
Ｂ４）ステップ２で求められた各画素領域の平均信号強度と、前記平均信号強度閾値とが比較されるステップと、
Ｂ５）前記領域の平均信号強度が、ステップ４の閾値より大きい場合、画素は追跡されるべきランドマークとして定義され、追跡されるべきランドマークリストに追加されるステップとから成る自動特徴選択ステップが実行されることを特徴とする方法。
【請求項２】
前記平均信号強度閾値が実験によって決定され、試験セット又はサンプル画像に対し前記方法を実行し、得られた異なる閾値の結果を比較することによって決定されることを特徴とする請求項１記載の方法。
【請求項３】
第２画像内の追跡されるべきランドマークに選択された画像の画素が、所謂ルーカス＆カナダ自動特徴追跡方法又はアルゴリズムによって自動的に追跡されることを特徴とする請求項１又は２記載の方法。
【請求項４】
請求項３記載の２次元画像にのみ導入されたルーカス＆カナダアルゴリズムの初期の実行が、２次元画像に代わって３次元画像にまで拡大され、
ａ）３次元画素により形成された第１の３次元画像ボリュームＩが提供されるステップと
、
ｂ）各目標画素に対して各画素を取り囲み、調整可能な例えば前記目標画素を中心に３×３×３アレイの画素を有する画素領域が定義されるステップと、
ｃ）デカルト座標系を定義し、前記デカルト座標系の各画素に対し相対的に、隣接画素に対して相対的に第１画像ボリュームの各画素の強度傾斜によって形成される下記式で定義される所謂傾斜ボリュームが計算されるステップと、
ｄ）各目標画素に対する以下の式で定義される所謂傾斜マトリックスが計算されるステップと、
ｅ）３次元画像ボリュームの各目標画素の各傾斜マトリックスに対し、以下の式を適用することで、最小固有値が計算されるステップと、
定義
傾斜マトリックスの最小固有値の計算
ｆ）最小固有値の閾値が定義されるステップと、
ｇ）対応する傾斜マトリックスが、
ｉ）固有値が閾値より大きい
ii）固有値が固有値のすべての最小値に対する最大値の割合より小さくないこと
iii）もし特徴として他の画素が、定義された選択された特徴としてまた維持された目標画素周りの画素領域内に存在する場合、傾斜マトリックスがより大きな最小固有値を有する画素のうちの一つのみが、特徴として選択され、他の一つは、追跡可能な特徴のリストから破棄される
という条件を満足する最小固有値を有する画像Ｉの追跡すべき特徴としての各画素が選択されるステップと、からなる生物医学的横断面画像ボリュームに適用されることを特徴とする方法。
【請求項５】
特徴として選択された画素の隣接領域の平均信号強度に対する閾値が約１０にセットされていることを特徴とする請求項４記載の方法。
【請求項６】
第１画像ボリュームＩの特徴として選択された各画素に対し、前記特徴が第２の遅い時間に同じ対象物の別の画像Ｊの画素アレイ内の位置に対して相対的に追跡され、前記追跡は、所謂ルーカス＆カナダの特徴追跡のピラミッド型実行によって実施されることを特徴とする請求項１ないし５のいずれか1記載の方法。
【請求項７】
画像Ｉの３次元画素アレイに対応する画像ボリュームＩ内のポイントとして、ポイント
ｕを、画像Ｊの3次元画素アレイに対応する画像ボリュームＪの対応するポイントをｖと定義し、これらのポイントが画像Ｉの特徴として選択された画素の座標を有するステップ
と、
当初のボリュームを表す基部（Ｌ＝０）と共に、レベルＬ＝０、．．．ｍまでのボリュームを表す、画像ボリュームＩとＪの２つの理想的なピラミッドＩLとＪLとが構築されるステップと、
直接隣接する領域の画素の一定量を一つの平均値と組み合わせることで、各追従する階のボリュームのサイズが減少されるステップと、
もっとも高い階Ｌmで初期化される所謂包括的ピラミッド型移動ベクトルｇを定義するステップと、
ｉ）ｕLとしてのポイントｕの位置を定義し、以下の式で前記位置を計算するステップ
と、ここでＬはピラミッドの実際の階又はレベルでの値であり、ここでｐは実際のボリューム位置であり、
ii）ボリューム傾斜のための以下の式に従って、デカルト座標ｘ、ｙ、ｚの各方向に対
してボリューム傾斜を計算するステップと、
iii）以下の式に従って傾斜マトリックスを計算するために傾斜ボリュームを用いるステップと、
ここでの値は、追跡された特徴を表現する目標画素に影響を与えるエリアサイズ又は画素の隣接領域を定義し、
iv）各レベルＬに対し、下の式で定義される反復ベクトルを初期化するステップと、
ｖ）kが１から最大カウント数Ｋまで、又は最小置換
まで、以下のボリューム差を計算するステップと、
vi）以下の式に従って、ミスマッチベクトルを計算するステップと、
vii）以下の式に従って、光学フローベクトルを定義するステップと、
ここでＧ-1はステップiii）で定義された逆傾斜マトリックスＧである。
viii）上で定義された光学フローベクトルを用いた双方向ベクトルを計算する以下のステップと、
これをレベルＬの最終光学フローと等しくする以下のステップと、
ix）前記ステップｉ）からviii）を各レベルで最終レベルＬ＝０になるまで繰り返し、下記式で定義される最終光学フローベクトルに達するステップと、
ｘ）追跡されるべき特徴を現すポイントｕに対応するボリュームＪのポイントｖの座標
を下の式によって定義するステップと、
xi）画像Ｉ内で選択された特徴に対応する各画素に対して、上記方法を繰り返すステ
ップと、
からなる所謂ルーカス＆カナダの特徴追跡のピラミッド型実行であることを特徴とする請求項６記載の方法。
【請求項８】
逆光学フローベクトルを画像Ｉの選択された特徴に対応する画素のポイントｕに対応する特徴を現す画素に対応して識別される画像Ｊの各ポイントｖに加えることによって、 又は光学フローベクトルをそれぞれ画像Ｉの選択された特徴として識別される画素に対応したポイントｕに加えることによって、２つの画像ＩとＪの登録が実行されることを特徴とする請求項７記載の方法。
【請求項９】
造影剤強化診断イメージングのための方法との組み合わせ、特に造影剤強化ＭＲＩ又は超音波イメージングを提供することを特徴とする方法であって、
ａ）ＭＲＩ又は超音波撮影で得られた、少なくとも第１と第２のデジタル画像又は同一対象物の断面画像セットを提供するステップであって、
前記画像が画素の2次元又は３次元アレイによって形成され、
同一組織又は同一組織帯又は同一の解剖学上の部分の走査が造影剤の存在下で第２の時間又はいずれかの以後の時間に実行され、
ｂ）第１画像又は画像セットの内部で、特定数のランドマークや所謂特徴を特定数の画素を選択することによって定義し、追跡すべき前記特徴のリストを生成するステップと、
ｃ）特徴として選択された各画素に対し、第１画像から第２画像まで又は画像セットから前記光学フローベクトルを定義することによって、第１画像から第２画像で特徴として選択された各画素又は画像セットの位置を追跡するステップと、
ｄ）前記逆光学フローを第２画像の画素又は画像セットに加えることで、第１又は第２画像又は画像セットを登録するステップと
からなる請求項１ないし８のいずれか１に記載の方法。
【請求項１０】
Ｂ１）第１画像又は第一断面画像セットの各画素周りの画素近傍領域を定義するステップと、
ここで前記画素近傍領域は、限定された数の画素で構成され、
２次元近傍領域が、単一画像の場合に選択され、
３次元近傍領域が、断面画像セットの場合に選択され、
Ｂ２）各画素に対し、前記画素の近傍領域の全画素の平均信号強度を定義するステップと、
Ｂ３）平均信号強度閾値を定義するステップと、
Ｂ４）ステップＢ２で求めた平均信号強度と、前記平均信号強度閾値を比較するステップと、
Ｂ５）前記近傍領域の平均信号強度がステップ４で求めた閾値より高い場合、画素が追跡されるべき特徴として定義され、追跡されるべき特徴のリストに追加されるステップと、
からなる更なる自動特徴選択ステップが、特徴選択ステップの後、特徴追跡ステップの前に実行されることを特徴とする請求項１ないし９のいずれか１に記載の方法。
【請求項１１】
２つのデジタル画像又は画像ボリュームを登録するためのコンピュータプログラムであって、請求項１ないし１０のいずれか１に記載の方法がコード化され、以下の擬似コードで特徴付けられるコンピュータプログラム。
【図１】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図２】
【図１１】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図２】
【図１１】
【公表番号】特表２００８−５３８０８０（Ｐ２００８−５３８０８０Ａ）
【公表日】平成２０年１０月９日（２００８．１０．９）
【国際特許分類】
【出願番号】特願２００７−５５３５９８（Ｐ２００７−５５３５９８）
【出願日】平成１８年１月３１日（２００６．１．３１）
【国際出願番号】ＰＣＴ／ＥＰ２００６／０５０５７１
【国際公開番号】ＷＯ２００６／０８２１９２
【国際公開日】平成１８年８月１０日（２００６．８．１０）
【出願人】（５０３４２７３５９）ブラッコ イメージング ソチエタ ペル アチオニ (19)
【氏名又は名称原語表記】ＢＲＡＣＣＯ ＩＭＡＧＩＮＧ Ｓ．Ｐ．Ａ．
【Ｆターム（参考）】
【公表日】平成２０年１０月９日（２００８．１０．９）
【国際特許分類】
【出願日】平成１８年１月３１日（２００６．１．３１）
【国際出願番号】ＰＣＴ／ＥＰ２００６／０５０５７１
【国際公開番号】ＷＯ２００６／０８２１９２
【国際公開日】平成１８年８月１０日（２００６．８．１０）
【出願人】（５０３４２７３５９）ブラッコ イメージング ソチエタ ペル アチオニ (19)
【氏名又は名称原語表記】ＢＲＡＣＣＯ ＩＭＡＧＩＮＧ Ｓ．Ｐ．Ａ．
【Ｆターム（参考）】
[ Back to top ]