(57)【特許請求の範囲】
【請求項１】
マスク像、および不透明化の前に得られた、前記マスク像に関連したそれぞれのデータおよび造影剤ボーラスの注入後に取得収集された造影像を使用するイメージング技術におけるコンピュータ支援の３次元画像識別照合方法において、
順次連続する軸方向ＣＴマスク像および造影像をそれぞれの等方性の３Ｄ体積においてリサンプリングし、
前記マスク体積における３Ｄ特徴点を選択し、
前記造影体積における対応関係量ないし対応点を求め、
生じる粗３Ｄ画像フローベクトルを逐次ランダムアルゴリズムによって処理しかつ対の少なくとも１つの前以て決められた百分率によって取り決められている動きパラメータ、すなわち変位および回転を最小２乗法で最適化して計算して、患者の動きが発見されるようにし、
患者の動きが発見された後に、前記マスク体積を相応に変換しかつ該マスク体積を前記造影体積からサブトラクションし、かつ
結果的に生じた体積を再構成しかつ表示する
ことを特徴とするコンピュータ支援の３次元画像識別照合方法。
【請求項２】
順次連続する軸方向ＣＴマスク像および造影像をリサンプリングする前記ステップは次のステップから成る、すなわち
スライス間ボクセルに対する３線補間を実施し、かつ
前記ボクセルに対して直交方向においてサブサンプリングを実施する
請求項１記載のコンピュータ支援の３次元画像識別照合方法。
【請求項３】
順次連続する軸方向ＣＴマスク像および造影像をリサンプリングする前記ステップが、画像識別照合アルゴリズムがＣＴＡアプリケーションにおいて所望される計算効率のレベルを実現することができるように画像データを低減する
請求項１記載のコンピュータ支援の３次元画像識別照合方法。
【請求項４】
前記画像識別照合アルゴリズムは、関心のある（注目）３Ｄ特徴点を突き止めるために、サブサンプリングを実施する前記ステップにおいてサブサンプリングされた関心のある点に対応付けられたサブサンプリングされたデータのみを使用する
請求項１記載のコンピュータ支援の３次元画像識別照合方法。
【請求項５】
前記関心のある点が突き止められた後、前記アルゴリズムは前記サブサンプリングされたデータを捨てる
請求項４記載のコンピュータ支援の３次元画像識別照合方法。
【請求項６】
前記サブサンプリングされたデータが捨てられた後、前記アルゴリズムは、前記最初に述べた、続いて行われる識別照合プロセスにおける原画像に戻って参照する
請求項５記載のコンピュータ支援の３次元画像識別照合方法。
【請求項７】
マスク像、および不透明化の前に得られた、前記マスク像に関連したそれぞれのデータおよび造影剤ボーラスの注入後に取得収集された造影像、および該造影像に関したそれぞれのデータを使用するイメージング技術におけるコンピュータ支援の３次元画像識別照合方法において、
順次連続する軸方向ＣＴマスク像および造影像をそれぞれの等方性の３Ｄ体積においてリサンプリングし、
前記マスク体積における３Ｄ特徴点を選択し、
前記造影体積における対応関係量ないし対応点を求め、
生じた粗３Ｄ画像フローベクトルを逐次ランダムアルゴリズムによって次のように処理する、すなわち
（ａ）所望の最大個別残留誤差εおよび所望のサイズｓによってスタートし、
（ｂ）Ｓを形成するために前記マスク体積の点の集合Ｖから小さな数の点、Ｎｌをランダムに選択しかつＳ′を形成するために造影体積の点の集合Ｗから相応の点を選択し、
（ｃ）単位４元数を使用して集合Ｓおよび集合Ｓ′に対する回転マトリクスＲおよび変位ベクトルｔを計算しかつ前記最大個別残留誤差がεより大きいならば、もはやこのようなケースでなくなるまで該ステップを繰り返し、
（ｄ）ＶおよびＷにおける残っている点から決まっている数の新しい点をランダムに取出し、かつ新しい変換パラメータを計算しかつ前記誤差制約条件が再び満たされたならば、これらの点をＳおよびＳ′に付加しかつそうでなければ、該ステップを繰り返し、
（ｅ）ＳおよびＳ′のサイズが、点がＶ′←ＳおよびＷ′←Ｓ′で終わる、ｓより大きいかまたはｓに等しくなるまで前記ステップ（ｄ）を繰り返し、またはステップ３の繰り返しの前以て決められた数、謂わばＴ１の後に、Ｓのサイズがｓに達していないならば、ステップ（ｂ）で再スタートし、または所定のεに対して、所望のサイズｓが前以て決められた数、Ｔ２の再スタート後に得られなかったならば、新しいεでステップ（ｂ）から再スタートする
ことを特徴とするコンピュータ支援の３次元画像識別照合方法。
【請求項８】
前記の回転マトリクスおよび変位ベクトルを計算した後、幾何学的変換を実施する
請求項７記載のコンピュータ支援の３次元画像識別照合方法。
【請求項９】
前記マスク体積を、新しい変換パラメータに従って変換しかつ表示目的のために等方性の体積に対してリサンプリングする
請求項８記載のコンピュータ支援の３次元画像識別照合方法。
【請求項１０】
前記造影体積も、同じ次元の等方性の体積に対してリサンプリングする
請求項９記載のコンピュータ支援の３次元画像識別照合方法。
【請求項１１】
サブボクセル強度値を、３線補間によって計算する
請求項１０記載のコンピュータ支援の３次元画像識別照合方法。
【請求項１２】
リサンプリングされた造影体積の、幾何学的に変換されかつリサンプリングされたマスク体積によるサブトラクションが行われて、サブトラクションされたＤＳＡ体積が形成される再構成ステップを含む、
請求項１１記載のコンピュータ支援の３次元画像識別照合方法。
【請求項１３】
前記の再構成ステップは、既知のタイプの再構成器を使用して実行される、
請求項１２記載のコンピュータ支援の３次元画像識別照合方法。
【請求項１４】
前記再構成ステップを実行し、これによりそれぞれ異なった観測角度における投影に相応する可変数の画像が出力されるようにする、
請求項１３記載のコンピュータ支援の３次元画像識別照合方法。
【請求項１５】
再構成される結果を、順次連続するフレームの迅速な表示を行うためのアニメータプログラムを使用して再生する
請求項１４記載のコンピュータ支援の３次元画像識別照合方法。
【請求項１６】
マスク像、および不透明化の前に得られた、前記マスク像に関連したそれぞれのデータおよび造影剤ボーラスの注入後に取得収集された造影像、および該造影像に関したそれぞれのデータを使用するイメージング技術におけるコンピュータ支援の３次元画像識別照合方法において、
順次連続する軸方向ＣＴマスク像および造影像をそれぞれの等方性の３Ｄ体積においてリサンプリングし、
前記マスク体積における３Ｄ特徴点のセットを選択し、
前記造影体積における点のセットとの対応関係量ないし対応点を求め、
結果的に得られた粗３Ｄ画像フローベクトルを逐次ランダムアルゴリズムによって次のように処理する、すなわち
（ａ）所望の最大個別残留誤差εおよび所望のサイズｓによってスタートし、
（ｂ）Ｓを形成するために前記マスク体積の点の集合Ｖから小さな数の点、Ｎｌをランダムに選択しかつＳ′を形成するために前記造影体積の点の集合Ｗから相応の点を選択し、
（ｃ）単位４元数を使用して集合Ｓおよび集合Ｓ′に対する回転マトリクスＲおよび変位ベクトルｔを計算しかつ前記最大個別残留誤差がεより大きいならば、もはやこのようなケースでなくなるまで該ステップを繰り返し、ここで単位４元数は単位回転を表し、
（ｄ）ＶおよびＷにおける残っている点から決まっている数の新しい点をランダムに取出し、かつ新しい変換パラメータを計算しかつ前記誤差拘束が再び満たされたならば、これらの点をＳおよびＳ′に付加しかつそうでなければ、該ステップを繰り返し、
（ｅ）ＳおよびＳ′のサイズが、点がＶ′←ＳおよびＷ′←Ｓ′で終わる、ｓより大きいかまたはｓに等しくなるまで前記ステップ（ｄ）を繰り返し、またはステップ３の繰り返しの前以て決められた数、謂わばＴ１の後に、Ｓのサイズがｓに達していないならば、ステップ（ｂ）で再スタートし、または所定のεに対して、所望のサイズｓが前以て決められた数、Ｔ２の再スタート後に得られなかったならば、新しいε（すなわちε←ε＋Δε）でステップ（ｂ）から再スタートする
ことを特徴とするコンピュータ支援の３次元画像識別照合方法。
【請求項１７】
ステップ（ｃ）において、前記マスク体積の、前記造影体積に対する変換のために回転マトリクスおよび変位ベクトルを得る
請求項１６記載のコンピュータ支援の３次元画像識別照合方法。
【請求項１８】
前記回転マトリクスおよび変位ベクトルを、誤差の２乗の和の最小化に従って計算する
請求項１７記載のコンピュータ支援の３次元画像識別照合方法。
【請求項１９】
前記回転マトリクスおよび変位ベクトルのそれぞれの中心点を得る
請求項１７記載のコンピュータ支援の３次元画像識別照合方法。
【請求項２０】
クロス（相互）−共分散マトリクスを計算し、該クロス−共分散マトリクスは、前記回転マトリクスおよび変位ベクトルにおける座標の積のそれぞれの和である要素を有している
請求項１９記載のコンピュータ支援の３次元画像識別照合方法。
【請求項２１】
前記単位４元数を、対称マトリクスの最も正の固有値に相応する単位固有ベクトルとして得る
請求項２０記載のコンピュータ支援の３次元画像識別照合方法。
【請求項２２】
前記対称マトリクスは４×４マトリクスであり、該マトリクスのそれぞれの要素は、前記回転マトリクスおよび変位ベクトルにおける座標の積のそれぞれの和である
請求項２１記載のコンピュータ支援の３次元画像識別照合方法。
【請求項２３】
それから最適な変位ベクトルを、前記造影体積における前記点のセットの前記中心点（Ａ）と前記マスク体積および前記回転マトリクスにおける３Ｄ特徴点の前記セットの前記中心点の積（Ｂ）との差として導出する
請求項２２記載のコンピュータ支援の３次元画像識別照合方法。
【請求項２４】
動きパラメータＲおよびｔを計算するステップ（ｃ）を次のステップに従って実施する：
【数１】
0003803404.tif 000017
４元数に基づいたアルゴリズムを使用し、ここにおいて単位４元数は、ｑ０≧０およびｑ０２＋ｑ１２＋ｑ２２＋ｑ３２＝１とする場合の１つの４ベクトルｑＲ＝［ｑ０ｑ１ｑ２ｑ３］ｔでありかつ単位回転を表すための方法であり、この場合Ｒを次のように計算し：
【数２】
0003803404.tif 000018
前記セットＶの中心点μｖを
【数３】
0003803404.tif 000019
に従って計算し、
３×３のクロス（相互）−共分散マトリクスＭを次のように計算し：
【数４】
0003803404.tif 000020
該マトリクスの要素は前記ＶおよびＷにおける座標の積の和であり、Ｍは次のように書き表すことができ：
【数５】
0003803404.tif 000021
最適な変位ベクトルを計算し：
【数６】
0003803404.tif 000022
ただし関心のある単位４元数は４×４対称マトリクスＮの最も正の固有値に相応する単位固有ベクトルであり、ただし
【数７】
0003803404.tif 000023
請求項７記載のコンピュータ支援の３次元画像識別照合方法。
【発明の詳細な説明】
【０００１】
【産業上の利用分野】
本発明は、ヘリカルスキャン方式のコンピュータ支援断層撮影血管造影法(スパイラル・コンピューテド・トモグラフィー・アンギオグラフィー)（computed tomography angiography＝ＣＴＡ）に関し、さらに特定すれば、骨像成分の抑圧およびそれに関連した、血管像成分の保持に対するＣＴＡ画像処理に関する。
【０００２】
【従来の技術】
ヘリカルスキャン方式のコンピュータトモグラフィ援用断層撮影血管造影法（ＣＴＡ）は、数秒内に大量のデータの収集を実現する、血管イメージングに対する新規な、所要穿刺が最小限ですむ（minimally invasive）技術である。M. A. F.ら‥E. H. Dillon, M. S. van Leeuwen,“Spiral ct Angiography” AJR, vol. 160, pp. 12731278, June 1993 参照。専門文献において報告されているＣＴＡに関するリサーチのうちの多くは、或る方法または別の方法における３次元の（３Ｄ）編集／分割技術を使用している。目的は、ＣＴ画像データにおける骨像に関するデータをできるだけ多く抑圧しかつ関心のある血管に関する画像データを含んでいるようにすることである。しかし、編集により、情報の損失および歪に対する可能性が生じかつしばしば平均で１５ないし３０分の操作時間を必要とする。また編集の品質は、主観的な経験にも相当依存している。
【０００３】
２次元の（２Ｄ）デジタルサブトラクション技術（一般にはデジタルサブトラクションアンギオグラフィー（Digital Subtraction Angiography＝ＤＳＡ）として知られている）はこの分野において従来より周知である。
【０００４】
本発明の１つの様相によれば、デジタルサブトラクション技術が３次元において利用される。
【０００５】
本発明の別の様相によれば、新規な３Ｄ画像識別照合（registration）アルゴリズムが３Ｄ ＤＳＡのために導入される。３Ｄ ＣＴ画像を離散的な軸方向位置においてサンプリングされた２Ｄ画像系列として取扱いかつ各スライスの動きアーチファクトを独立に補正する既存の技術とは異なって、本発明によれば、患者の大域的な動きが局所的な無意識の／滑らかな（柔軟な）動きに結び付いているときですら、それはサブトラクションの前に補正される。実験結果によれば、数セットの、臨床的なヘリカルスキャン方式によるＣＴデータに基づいた本発明によるこのアルゴリズムの使用の有効性が確認されている。
【０００６】
上述したように、ヘリカルスキャン方式ＣＴ血管造影法（ＣＴＡ）は、数秒内に大量のデータを収集することができるようにする、血管撮影に対する新規な、所要穿刺が最小限ですむ技術である。１  このデータ収集速度により、しばしば１回の息こらえにおいて、患者の少ない動きで実施することができる研究が可能になる。静脈の造影剤の静脈内（ＩＶ）ボーラスが注入されかつ適正に時間を置かれるならば、ヘリカルＣＴスキャンは、患者の解剖学的構造の大きな体積における動脈相を捕捉することができ、これにより血管の管腔、狭窄および病変の優れた視覚化が可能になる。それからこれらのデータは、３Ｄ視覚化(イメージング)技術（例えば体積再構成（＝volume rendering）、最大強度投影（maximum intensity projection＝ＭＩＰ）、およびシェーディングサーフェス表示）を使用して表示し、従来の血管造影撮影法の画像に類似した動脈の解剖学的構造の画像を得ることができる。
【０００７】
専門文献において報告されているＣＴＡにおける多くのリサーチは、或る方法または別の方法における編集技術を使用している。例えば、J. E. J. et al‥M. P. Marks, S. Napel,“Diagnosis of carotid artery disease: Preliminary experience with maximu-intensity-projection spiral ct angiography,” AJR, vol. 160, pp. 1267-1271, june 1993 、G. D. R. et al‥S. Napel, M. P. Marks,“Ct angiography with spiral ct and maximum intensity projection,”Radiology, vol. 185, pp. 607-610, Nov.1992 、D. M. C. et al., R. B. Schwartz, K. M. Jones,“Commomcarotid artery bifurcation: Evaluation with apiral ct,”Radiology, vol.185, no.2, pp.513-519, 1992 が参考になる。
【０００８】
上述したように、この目的は、ＣＴデータにおける骨像データをできるだけ多く抑圧しかつ関心のある血管に関する画像データを含んでいるようにすることである。骨およびカルシウム沈積物に関する画像データの除去のために、ＭＩＰ再構成法（reconstruction）は、血管の管腔内の造影剤コラムを一層明瞭に表示することができる。しかしこの編集手法は幾つかの制限を有している。第１に、編集により、情報の損失または歪に対する可能性が生じる。このことは、連結性アルゴリズムが、血管と物理的に接触している骨を排除するために著しく低いしきい値を使用しているときに起こる。このことは、編集プロシージャが、血管の管腔に接触する血管壁カルシウム沈積を抑圧しようとするときも起こり得る。第２に、編集には典型的には、平均１５ないし３０分の操作時間が必要でありかつ、場合によっては何時間もかかることすらある。さらに、多くの場合編集の品質は非常に主観的でありかつ個人の経験に負うところが大きい。例えば、壁内のカルシウムを管腔内の造影剤から区別できるようにするために、強調されていない画像に対する基準が必要である。
【０００９】
ヘリカル方式ＣＴＡに対する択一的な手法は、従来のデジタルサブトラクション技術（一般にＤＡＳとして知られている）を２次元から３次元に拡張することである。ＤＡＳでは２つの画像セットが使用される：データの第１セットは、不透明化の前に得られるマスク像であり、一方造影像と称されるデータの第２セットは、造影剤ボーラスの注入後に収集される。マスクに関連した画像が、造影に関連した画像からサブトラクションされる。マスク像および造影像における血管の位置が一致すれば、サブトラクションにより純然たる樹木状の血管像が形成される。しかしこの理想的な状況は実際の臨床では実現するのが常に難しい。事実、数多くの研究には、患者の動きがこのサブトラクション技術を大きく制限しておりかつ結果として生じるアーチファクトを除去するために種々の方法が示唆されていることが示されている。データ収集の間、患者の動きを制限するために、頭部保持装置を使用することができ、一方後処理の期間に患者の動きを補正するために、時間フィルタリングおよびより一般的には画像識別照合が適用される。例えば次の文献を参照されたい。B. B. J., L. Erikson, T. Greitz, T. Ribbe, and L. Widen,“Head fixation device for reproducible position alignment in transmission ct and pet,”J. Comp. Assis. Tomgr., pp.136-141, 1981。
【００１０】
既存のＤＡＳ技術は３Ｄ ＣＴ画像を、離散的な軸線方向の位置においてサンプリングされた２Ｄ画像として取り扱いかつ動きアーチファクトはマスク像および造影像の各スライス間の識別照合処理を通して独立して補正される。例えば、J. J. G. rt al‥J. M. Fitzpatrick, D. R. Pickens III,“Technique for automatic motion correction in digital subtraction angiography”Optical Engineering, vol. 26, pp. 1085-1093, Nov. 1987 が参考になる。しかしここにおいて、患者の動きは普通３次元的でありかつそれは３次元の動きパラメータを発見することによって補正されるべきであることが認められている。
【００１１】
大抵の従来の３Ｄ画像識別照合技術は、ＣＴ，ＰＥＴ（Positron Emission Tomography），ＳＰＥＣＴ（Single Positron Emission Tomography），ＭＲＩ（Magnetic Resonance Imaging）および超音波イメージング技術のような種々の形式から得られる３Ｄ画像の融合において使用される。知られている限りでは、これら技術のいずれも、３Ｄ ＣＴＡの問題に適用されていない。
【００１２】
伝統的には、外部マーカ、定位的なフレームおよび解剖学的な標識が、画像識別照合を必要とする用途に使用されている。例えば、G. Chen, M. Kessleer, and S. Pitluck,“Structure transfer in the 3d medical imaging studies,”in Computer graphics,(TX), pp. 171-175, National Computer Graphics Association, 1985 が参考になる。
【００１３】
これらの方法は、コントロールされた環境で実施されなければならずかつ人体の解剖学的な特徴認識についての専門知識を要求する。脳像に対する半自動的な空間識別照合スキーマが、Bartoo および Hanson によって開発され、ここでは基準点として解剖学的な対象物の中心（セントロイド）が使用される。参照［８］G. Bartoo and W. Hanson“Multi-modality image registration using centroid mapping,”IEEE Engineering in Medicine and Biology Society 11th Annual International Conference, 1989。これらの技術では、主要な解剖学的な対象物の中心が、歪の存在する場合に近似的に同じ中心を維持しかつ手動の識別照合の結果に匹敵する結果をもたらすということが前提になっている。
【００１４】
Pelizzari et al. は、マルチプルな脳像の３Ｄ識別照合に対するサーフェスマッチング技術を開発した。ここでは、頭部の外面が、それは脳のすべてのスキャンにおいて可視状態であるので、基準面として使用される。C. Pelizzari, G. Chen, D. Spelbring, R. Weichselbaum, snd C. Chen,“Accurate three-dimentional registration of ct, pet, and/or mr images of the brain,”Journal of Computer Assisted Tomography, vol. 13, pp. 0-26, January/February 1989 参照。この基準面により、患者固有の座標系が確定される。この３Ｄ識別照合は、それぞれのスキャンのシリアルスライスに基づく種々異なったアウトラインの輪郭間の最良のマッチングを最小２乗の手法において逐次近似的に（反復的）に見つけ出すことによって行われる。それから動きパラメータ、すなわち回転、変位およびスケーリング係数が計算される。それらの方法は、頭部モデル上に位置する部分面を最小２乗適合プロシージャから排除することによるこれら部分面の或る程度の消失を許容することができるが、逐次近似特性は、迅速な収束に対する申し分ない初期推定を必要とする。“ステアリング”と呼ばれる人間の介入は、適正な収束速度を保証するために変換パラメータを調整することが要求される。参考になるのは、C. Pelizzari, G. Chen, D. Spelbring, R. Weichselbaum, snd C. Chen,“Accurate three-dimentional registration of ct, pet, and/or mr images of the brain,”Journal of Computer Assisted Tomography, vol. 13, pp. 0-26, January/February 1989。
【００１５】
Alpert et al は、２セットの３Ｄ画像データを識別照合するために主軸変換技術を使用しており、ここではセットそれぞれが全体の脳体積をカバーしている。N. Alpert, J. Bradsshaw, D. Kennedy, and J. Correia,“The principal axes transformation-a method for image registration,”J Nucl Med, vol. 31, pp. 1717-1722, October 1990 参照。この手法では、それぞれのデータセットに対してマス（塊）の中心およびその主軸が計算される。マスの中心および３つの主軸は一緒に、１つの基準フレームを定める。それから、１つの基準フレームをどのように別のフレームに変換することができるかを見つけ出すことによって識別照合を行い得る。このアプローチでは、関心のある体積の手動のアウトライン化が必要であり、これに基づいてマスの中心および主軸の計算が行われる。この手法では計算の簡素化およびスピード化が実現されるが、頭部全体が両方の体積において存在していなければならないので、厳しい制限を受けている。上に挙げたN. Alpert等による刊行物を参照されたい。
【００１６】
また、Moshfeghi and Rusinek により、類似の３Ｄ識別照合法が提供されている。ここでは、固有値解析に対して、主軸の位置決めのためにまばらなマトリクスが使用されかつそれから回転、変位およびスケーリングパラメータが計算される。ここでは画像中に局所的な歪は現われないものと仮定している。M. Moshfeghi and H. Rusinek,“Three-dimensional registration of mulltimodality medical images using the principal axes technique,”Philips Journal Research, vol. 47, pp. 81-97, 1992。
【００１７】
最近、Thirion et al. により、３Ｄ画像の両セットから抽出される最高線（クレストライン）のマッチングに基づいた自動３Ｄ識別照合アルゴリズムが提供されている。最高線は、３Ｄ対象物の表面における高い曲率点に対応するラインである。この方法はまず、関心のある対象物面を抽出するために３Ｄ画像を分割しかつそれから最高線を自動的に抽出するためにマーチングラインを使用する。２つのデータセットからの最高線は、幾何学的な切断に基づいてマッチングされる。A. Gueziec and N. Ayache,“Smoothing and matching of 3d-space curves,”in Proceedings of the Second European Conference on Computer Vision, May 1992 参照。このプロセスの終わりで、非常に多数の点をマッチングする変換が行われる。S. B. J.-P. Thirion, O. Monga and etc.,“Automatic registration of 3d images using surface curvature,”in SPIE Mathematical Methods in Medical Imaging, vol. 1768, pp. 206-216, 1992 参照。
【００１８】
【発明の開示】
本発明は、３Ｄ ＣＴ ＤＳＡにおけるこのような問題の解決に適用されるアルゴリズムを使用した３次元画像識別照合方法を開示する。本発明は、識別照合における不完全なまたは部分的な体積（ボリューム）データを扱いかつ患者の大域的な動きが局所的に無意識な／柔軟な動きと結び付いているときですら、サブトラクションの前に患者の大域的な動きを補正することができる。実験結果は、数個の臨床ヘリカルＣＴデータに基づいたこのアルゴリズムの効能を示している。
【００１９】
本発明の様相によれば、マスク像、および不透明化の前に得られた、前記マスク像に関連したそれぞれのデータおよび造影剤ボーラスの注入後に取得収集された造影像を使用するイメージング技術におけるコンピュータ支援の３次元画像識別照合方法は、次のステップから成っている：順次連続する軸方向ＣＴマスク像および造影像をそれぞれの等方性の３Ｄ体積においてリサンプリングし；前記マスク体積における３Ｄ特徴点を選択し；前記造影体積における対応関係量ないし対応点を求め；生じた粗３Ｄ画像フローベクトルを逐次ランダムアルゴリズムによって処理しかつ対の少なくとも１つの前以て決められた百分率によって取り決められている最良の動きパラメータ、すなわち変位および回転を最小２乗法において計算して、患者の動きが発見されるようにし；患者の動きが発見された後に、前記マスク体積を相応に変換しかつ該マスク体積を前記造影体積からサブトラクションし；かつ結果的に生じた体積を再構成しかつ表示する。
【００２０】
本発明の別の様相によれば、順次連続する軸方向ＣＴマスク像および造影像をリサンプリングするステップは、スライス間ボクセルに対する３線補間を実施し、かつボクセルに対して直交方向、すなわちＸおよびＹ方向においてサブサンプリングを実施する。
【００２１】
本発明の別の様相において、順次連続する軸方向ＣＴマスク像および造影像をリサンプリングするステップは、画像識別照合アルゴリズムがＣＴＡアプリケーションにおいて所望される計算効率のレベルを実現することができるように画像データを低減する。
【００２２】
本発明のさらなる様相において、画像識別照合アルゴリズムは、関心のある３Ｄ特徴点を突き止めるために、サブサンプリングを実施する前記ステップにおいてサブサンプリングされた関心のある点に対応付けられたサブサンプリングされたデータのみを使用する。
【００２３】
本発明のさらなる様相によれば、関心のある点が突き止められた後、アルゴリズムは前記サブサンプリングされたデータを捨てる。
【００２４】
本発明のさらなる様相によれば、サブサンプリングされたデータが捨てられた後、アルゴリズムは、最初に述べた、続いて行われる識別照合プロセスにおける原画像に戻って参照する。
【００２５】
本発明のさらに別の様相において、マスク像、および不透明化の前に得られた、前記マスク像に関連したそれぞれのデータおよび造影剤ボーラスの注入後に取得収集された造影像、および該造影像に関したそれぞれのデータを使用するイメージング技術におけるコンピュータ支援の３次元画像識別照合方法は次のステップを有している：順次連続する軸方向ＣＴマスク像および造影像をそれぞれの等方性の３Ｄ体積においてリサンプリングし；マスク体積における３Ｄ特徴点を選択し；造影体積における対応関係量ナイシ対応点を求め；生じた粗３Ｄ画像フローベクトルを逐次ランダムアルゴリズムによって次のように処理する、すなわち（ａ）所望の最大個別残留誤差εおよび所望のサイズｓによってスタートし、（ｂ）Ｓを形成するためにＶから小さな数の点、謂わばＮｌをランダムに選択しかつＳ′を形成するためにＷから相応の点を選択し、（ｃ）単位４元数を使用してＳおよびＳ′に対する回転マトリクスＲおよび変位ベクトルｔを計算しかつ前記最大個別残留誤差がεより大きいならば、もはやこのようなケースでなくなるまで該ステップを繰り返し、（ｄ）ＶおよびＷにおける残っている点から決まっている数の新しい点をランダムに取出し、かつ新しい変換パラメータを計算しかつ誤差制約条件が再び満たされたならば、これらの点をＳおよびＳ′に付加しかつそうでなければ、該ステップを別のやり方で行い、（ｅ）ＳおよびＳ′のサイズが、点がＶ′←ＳおよびＷ′←Ｓ′で終わる、ｓより大きいかまたはｓに等しくなるまで前記ステップ（ｄ）を繰り返し、またはステップ３の繰り返しの前以て決められた数、謂わばＴ１の後に、Ｓのサイズがｓに達していないならば、ステップ（ｂ）で再スタートし、または所定のεに対して、所望のサイズｓが前以て決められた数、Ｔ２の再スタート後に得られなかったならば、新しいε（すなわちε←ε＋Δε）でステップ（ｂ）から再スタートする。
【００２６】
本発明のさらに別の様相において、回転マトリクスおよび変位ベクトルを計算した後、幾何学的変換を実施する。
【００２７】
本発明のさらになお別の様相において、マスク体積を、新しい変換パラメータに従って変換しかつ表示目的のために等方性の体積に対してリサンプリングし、造影体積も、同じ次元の等方性の体積に対してリサンプリングし、サブボクセル強度値を、３線補間によって計算し、リサンプリングされた造影体積の、幾何学的に変換されかつサンプリングされたマスク体積によるサブトラクションによって、サブトラクションされたＤＳＡ体積を形成し、サブトラクションされたＤＳＡ体積を、市販のアプリケーション、the Ney/Fishman's volume renderer Iprender を使用して再構成し、再構成アルゴリズムにより、それぞれが異なった観測角度における投影に相応する可変数の画像を出力しかつ再構成された結果を、ｘ線技師が、そうしなければスタチックな観測からは使用可能でない、比較的多量の奥行きデータを合成することができるようにするために、順次連続するフレームの迅速な表示を行うためのアニメータプログラムを使用して再生（review）する。
【００２８】
本発明の別の様相において、方法は次のステップから成る：すなわちマスク体積における３Ｄ特徴点のセットを選択し；かつ造影体積における点のセットとの対応関係量ないし対応点を求める。
【００２９】
本発明のさらに別の様相において、マスク体積の、造影体積に対する変換のために回転マトリクスおよび変位ベクトルを得、回転マトリクスおよび変位ベクトルを、誤差の２乗の和の最小化に従って計算し、データ点の２つのセットの中心点を得、クロス−共分散マトリクスを計算し、該クロス−共分散マトリクスは、２つのデータセットにおける座標の積のそれぞれの和である要素を有しており、単位４元数を、対称マトリクスの最も正の固有値に相応する単位固有ベクトルとして得、対称マトリクスは４×４マトリクスであり、該マトリクスのそれぞれの要素は、２つのデータセットにおける座標の積のそれぞれの和であり、かつそれから最適な変位ベクトルを、造影体積における点のセットの中心点（Ａ）とマスク体積および回転マトリクスにおける３Ｄ特徴点のセットの中心点の積（Ｂ）との差として導出する。
【００３０】
【実施例】
次に本発明を図示の実施例につき図面を用いて詳細に説明する。
【００３１】
図１には、本発明の種々のステップが示されている。連続的な軸方向ＣＴ画像はまず、等方性の３Ｄ体積、すなわち等方性の画像解像力を呈する体積においてリサンプリングされる。それからマスク体積における３Ｄ特徴点が選択されかつ造影体積において対応関係量ないし対応点が求められる。それから結果として生じる粗３Ｄ画像フローベクトルが、対の少なくとも前以て決められた百分率において取り決められている最小２乗の手法において最良の動きパラメータ（変位および回転）を計算するために逐次近似ランダムアルゴリズムによって処理される。患者の動きが発見された後、マスク体積が相応に変換されかつ造影体積からサブストラクションされる。それから結果的にサブストラクションされた体積が再構成されかつディスプレイ上にて観測される。次にこの方法をより詳しく説明する。
【００３２】
本発明による第１のステップは、原マスク像を等方性データセットにリサンプリングすることである。これに関連して、１ボクセルは１つの３次元の画像単位成分、例えば単位立方体であり、２次元におけるピクセル、例えば単位正方形と類似している。スライス間ボクセルのために３線補間が使用され、一方ボクセルに対してｘ方向およびｙ方向においてサブサンプリングが行われる。このサンプリングプロセスにより画像データは低減されかつそれ故にＣＴＡアプリケーションにおいて所望されるレベルの計算効率を実現するために本発明によるアルゴリズムが許容される。
【００３３】
しかし、サブサンプリング技術を使用する従来の形式の手法とは異なって、本発明では、画像記述の成功に対して問題である精度のレベルを維持している。本発明による画像記述アルゴリズムでは、サブサンプリングされた点は、３Ｄ特徴点の場所を突き止めるためにのみ使用される。これらの関心ある点を突き止めた後、アルゴリズムは、サブサンプリングされたデータを捨てかつ常に、いずれの順次連続する記述プロセスにおいても原画像に戻って関連付ける。この様相もしくは態様は続く説明において一層明らかになる。
【００３４】
異なった３Ｄ特徴は、独自のマッチングを持つのが尤もらしい傾向にありかつそれ故に画像識別照合の精度を改善することになる。従来技術に報告されている大抵の特徴検出作業は、２次元画像に対してである。３次元の特徴は、視覚化するのが困難でありかつ場所を突き止めるために相当の時間を必要とする。Moravec によって表された、２Ｄから３Ｄへの関心のある演算子がここに展開されかつ迅速であるばかりでなく、非常に申し分なく働くことが認められている。H. P. Moravec,“Rover visual obstacle avoidance,”in Proceedings of the 7th International Conference on Artificial Intelligence, (Vancouver, B. C. , Canada), pp. 785-790. August 1981 参照。２次元の場合、関心のある演算子は、４つの可能な方向に対する各ピクセル位置における強度差の２乗の和（分散の推定）を計算する。関心のある値は、すべての分散値の和の最小値として定義されている。この演算子が３次元に拡張されるとき、それはそれぞれのボクセル位置において、図２に示されているように可能な１３の方向すべてにおける分散値の最小値であると定義されている、関心のあるインデックスを計算する。
【００３５】
この関心のあるインデックスｉは、数学的に、次のように表すことができる：
【００３６】
【数８】
0003803404.tif 000002
【００３７】
それからこのアルゴリズムは、最も高い関心のあるインデックスを有する前以て決められた数の点を特徴点とすべく選択する。全体の体積において一層均一に分配された特徴点を保証するために、局所最大値フィルタが、小さな体積内の特徴点のクラスタ（塊）化を回避するために使用される。図３には、異なった頭部分の種々の観測からの３つの２Ｄ画像スライスにおいて選択されたいくつかの特徴点が示されている、即ち：（左）軸方向像；（真中）冠状断層像（コロナル像）；（右）矢状断層像（サジタル像）。
【００３８】
【数９】
0003803404.tif 000003
【００３９】
【数１０】
0003803404.tif 000004
【００４０】
【数１１】
0003803404.tif 000005
【００４１】
【数１２】
0003803404.tif 000006
【００４２】
画像フローフィールドが計算された後、アルゴリズムは患者の動きパラメータ、すなわち変位および回転パラメータを推定しようとする。その場合動き推定は、逐次ランダムアルゴリズムを使用して実施される。臨床現場において、プロセス全体に１．５ないし２分のオーダがかかるので、患者が安定状態を維持することは難しい。食道の動きのような筋肉の収縮による局所的な無意識の動きが確かにある。
【００４３】
本発明によれば、画像フローベクトルの少なくとも１つの前以て決められた百分率によって取り決められている最良の動きパラメータを最小２乗の手法において計算する逐次ランダムアルゴリズムが使用されている。このパラメータ推定アルゴリズムは、R. M. Karp,“An introduction to randomized algorithms,” Discrete Applied Mathematic, vol. 34, pp. 176-210, 1991 に記載されているような、４元数手法およびおよびランダムアルゴリズムストラテジーに基づいている。B. K. P. Horn,“Closed-form solution of absolute orientation using unit quaterions,”Opt. Soc. Am., vol. 4, pp. 629-642, April 1987 参照。便利なように、４元手法の詳細が準備した参考文献の補遺に示してある。
【００４４】
【数１３】
0003803404.tif 000007
【００４５】
換言すれば、当該セットにおけるすべての点の残留誤差の最大値が最小化されるような点の可能な最大のセット。
【００４６】
この問題は、徹底的な探索を通して決定的に解決することができる。しかし、実施時間および記憶要求は極めて高いものになる。最適な解決法を推定する逐次ランダムアルゴリズムが提案される。このアルゴリズムは小さな残留誤差パラメータεによってスタートし、かつそれからεより小さい最大の個々の残留誤差でＶ′およびＷ′を見付けるようにする。その場合十分小さなεに対する解決法は、最適な解決法を近似することである。このアルゴリズムは次の通りである：
１．所望の最大個別残留誤差εおよび所望のサイズｓを以ってスタートする。
【００４７】
【数１４】
0003803404.tif 000008
【００４８】
３．ＶおよびＷにおける残留点から決められた数の新しい点をランダムに取上げ、かつ新しい変換パラメータを計算する。誤差制約条件がここでも満たされているならば、これらの点をＳおよびＳ′に付加する。その他の場合はこのステップを繰り返す。
【００４９】
４．（ａ）ＳおよびＳ′のサイズが、そこでは点がＶ′←ＳおよびＷ′←Ｓ′で終了するｓより大きいかまたはｓに等しくなるまで、ステップ３を繰り返し、または
（ｂ）ステップ３の繰り返しの何回か（Ｔ１）の後、Ｓのサイズがｓに達しないならば、ステップ２から再スタートし、または
（ｃ）所定のεに対して、所望のサイズｓが多数回の再スタート（Ｔ２）の後得られなかったならば、新しいε（即ちε←ε＋Δε）で以ってステップ２から再スタートする。
【００５０】
回転マトリクスおよび変位ベクトルが計算された後、幾何学的な変換が実施される。マスク体積はこれらパラメータに応じて変換されかつ同時に、表示目的で１つの等方性の体積に対してリサンプリングされる。造影体積も、同じ次元の等方性体積に対してリサンプリングされる。サブボクセル強度値は、３線補間によって計算される。リサンプリングされた造影体積の、幾何学的に変換され（かつサンプリングされた）マスク体積によるサブトラクションによって、サブトラクションされたＤＳＡ体積が形成される。
【００５１】
サブトラクションされたＤＳＡ体積は、市販のアプリケーション、the Ney/Fishman's volume renderer IPRENDER を使用して再構成される。このような体積の再構成技術は最初、Drebin et al. によって開発されかつ後にＣＴデータに応用された。R. A. Drebin, L. Carpenter, and P. Hanrahan,“Volume redering,”Computer Graphics (Proc. SIGGRAPH), vol. 22, pp. 65-74, August 1988 および D. R. Ney, E. K. Fishman, and D. Magid,“Volumetric rendering of computed tomography data: Principles and techniques,” IEEE Computer Graphics and Applications, pp. 24-32, March 1990 参照。
【００５２】
この再構成アルゴリズムは、それぞれが種々異なった観測角度における投影に相応する、可変数の画像を出力する。再構成された結果は、アニメータプログラムを使用して観測される。
【００５３】
順次連続するフレームの迅速な表示により、ｘ線技術者は、その他の場合には静的な観測から使用可能ではない、比較的多量の奥行き情報を合成することができる。
【００５４】
動き推定に対する４元数解決法は次の通りである。
【００５５】
【数１５】
0003803404.tif 000009
【００５６】
Ｒおよびｔは、誤差の２乗の和
【００５７】
【数１６】
0003803404.tif 000010
【００５８】
の最小化を前提として計算することができる。
【００５９】
動きパラメータを計算するために、４元数に基づいたアルゴリズムが使用される。このアルゴリズムがここで説明される。もっと詳細には、B．K．P．Horn，“Closed-form solutuon of absolute orientation using unit quaternions”Opt．Soc．Am．，vol．4，pp．629-642，April  1987  参照。
【００６０】
単位４元数は４つのベクトルｑＲ＝［ｑ０ｑ１ｑ２ｑ３］ｔであり、ただしｑ０≧０およびｑ０+２ｑ１+２ｑ２+２ｑ３+２＝１。単位回転数を表すための択一的な方法は次の通りである：
【００６１】
【数１７】
0003803404.tif 000011
【００６２】
【数１８】
0003803404.tif 000012
【００６３】
【数１９】
0003803404.tif 000013
【００６４】
が計算される。その要素は、ＶおよびＷにおける座標の積の和である。このマトリクスは、回転に対する最小２乗問題を解決するために要求される情報をすべて含んでいる。Ｍが次の式
【００６５】
【数２０】
0003803404.tif 000014
【００６６】
それから関心のある単位４元数は、４×４の対称マトリクスＮの最大の正の固有値に相応する単位固有ベクトルであり、ただし
【００６７】
【数２１】
0003803404.tif 000015
【００６８】
【数２２】
0003803404.tif 000016
【００６９】
本発明による実験的な使用例において、患者の頭部の比較的低い部分（上頸部から目の下の領域まで）の画像の２セットが、Siemens Somatom Plus Spiral Volumetric Computed Tomography(CT) scanners を使用して収集された。これは試し的な実行でしかなかったので、データ収集プロセスにおいて特別な注意は払われなかった。２つのスキャン間の待ち時間は、数分のオーダにあることが知られている。
【００７０】
各画像スライスのピクセルサイズはｘ方向およびｙ方向ともに０．２９６９７mmであり、３mmのスライス間分離を有する。両方の画像セットに対する解像力は、５１２ピクセル×５１２ピクセルでありかつ画像セット当たり全部で５７のスライスである。収集されたデータは０および２０４７間の強度値範囲を有している。患者は、頸部からウィリスの輪（大脳動脈輪）までスキャンされた。スキャンの低い方の半部は、多くの柔軟な動きが集中している頸部を含んでいる。
【００７１】
画像セットは、口内の２つの金歯のために５７のスライスのうち１２（スライス番号２０ないし３２）に重大な欠陥を有していることがわかっている。図４には、いくつかの２Ｄマスクおよび収集された造影像スライスが示されている。
【００７２】
本発明による手法では、原データの入力、サンプリングプロセス、３Ｄ画像識別照合、幾何学的な変換およびスブトラクションされた体積の最終的な出力を含んでいる、SUN SPARCstation １０における１２８×１２８×１４４の解像力でのＣＰＵ時間の１３分より短い時間がかかった。２５６×２５６×２８８の解像力では、一層長いＣＰＵ時間がかかり、その内大抵は、Ｉ／Ｏ、幾何学的な変換およびサンプリングに費やされる。特徴抽出、マッチングおよびパラメータ推定部分には、計算の煩雑さが主に、今日では両方の場合とも３００である使用される特徴点の数に依存しているので、大体同じＣＰＵ時間を要する。
【００７３】
再構成されサブトラクションされた体積は、アニメーションソフトウェアを使用してスクリーン上に表示される。アニメーションディスプレイのスナップショットは、スクリーンにて捕捉されかつ図５に示されているものである。同じ体積の、異なった角度からの別の観測が図７に示されている。比較として、同じ観測からのサブトラクションされた体積が図６および図８に示されている。これらは、サンプリングされたマスク体積をサンプリングされた造影体積から直接（識別照合なしに）サブトラクションすることによって得られたものである。
【００７４】
さらに、識別照合による改善を明らかにするために、図９には、サブトラクションされた体積の数個の２次元軸線方向スライスが示されている。強度値が暗ければ暗いほど、マスクおよび造影内の血管の位置はより良好に一致する。
【００７５】
識別照合された画像セットのＣＴＡ体積は、識別照合されないＣＴＡ体積に関する３Ｄ空間における大脳血管の視覚化を著しく改善した。にも拘らず開発されたＣＴＡ手法は一般に、スキャンされた頭部領域の大抵の部分における血管の高度な視覚化を提供した。
【００７６】
結果は明らかに、患者の動きの大多数が補正されたことを示している。アルゴリズムは、２つの植え込まれた金歯のためのイメージングアーチファクトの影響を回避することができる。しかし、（患者のえん下動作による）食道の筋肉収縮および（えん下動作に伴う可能性がある）下顎の動きが残る。
【００７７】
本発明による３Ｄ ＣＴＡに対するアルゴリズムにより、少なくとも次の利点が得られる：
Ａ．このアルゴリズムは、操作員支援による編集に対して要求される時間に殆ど匹敵する。
【００７８】
Ｂ．この識別照合アルゴリズムは、血管が骨に物理的に接触しているときですら、骨像を除去しかつ関心のある血管像を維持するのに成功する。
【００７９】
Ｃ．このアルゴリズムは、スキャニングプロセスの間の３Ｄの患者の動き並びにテーブルの動きの不正確さを処理することができる。
【００８０】
Ｄ．３Ｄの関心のある演算子は迅速でありかつ３Ｄのコーナ点の検出において非常に良好に作業する。
【００８１】
Ｅ．この識別照合アルゴリズムはその計算効率を改善するために原画像をサブサンプリングするにも拘らず、常に、計算の各段階において原画像に戻って関連付けて、所望のサブピクセル精度が得られるようにする。
【００８２】
Ｆ．大域的な患者の動きが常に、局所的な無意識／柔軟な動きに結合している一方、逐次近似ランダムアルゴリズムは、サブトラクションプロセスの前に患者の動きの大多数を補正することができる。
【００８３】
Ｇ．このアルゴリズムは、３Ｄ画像識別照合を要求するマルチモードデータ統合のような他の医療応用に容易に適合することができる。
【００８４】
本発明は、コンピュータ具体化によって実現するのに向いている。本発明を実施例を用いて説明してきたが、当業者に当然と思われる種々の変形および修正は、本発明の精神を逸脱しない限り実現可能である。この種の変形は、請求項によって定められている本発明の範囲内において行われている。
【図面の簡単な説明】
【図１】本発明の方法の外観をフローチャートにより示す図である。
【図２】本発明の理解に役立ちかつ３Ｄ特徴点を見つけ出すための１３の方向を表すベクトル表示を示す図である。
【図３】選択されたいくつかの特徴点を示す２Ｄ画像の図である。
【図４】スキャンされた体積の典型的な原マスク像および造影像スライスを示す図であり、（ａ）スライス５７（ｂ）スライス４９（ｃ）スライス４０（ｄ）スライス２７（２つの金歯によって惹起されているアーチファクトに注意）（ｅ）スライス１７（ｆ）スライス９（ｇ）スライス１。
【図５】３Ｄ画像識別照合およびデジタルサブトラクション後の再構成されたＣＴＡ体積の観測を示す図である。
【図６】３Ｄ画像識別照合なしの再構成されたＣＴＡデジタルサブトラクション体積の観測を示す図である。
【図７】３Ｄ画像識別照合およびデジタルサブトラクション後の再構成されたＣＴＡ体積の別の観測を示す図である。
【図８】３Ｄ画像識別照合なしの再構成されたＣＴＡデジタルサブトラクション体積の別の観測を示す図である。
【図９】いくつかの典型的な２Ｄ画像スライスに対する識別照合されない対識別照合された絶対デジタルサブトラクション差を示す図である。
【図１】
0003803404.tif 000024
【図２】
0003803404.tif 000025
【図３】
0003803404.tif 000026
【図４】
0003803404.tif 000027
【図５】
0003803404.tif 000028
【図６】
0003803404.tif 000029
【図７】
0003803404.tif 000030
【図８】
0003803404.tif 000031
【図９】
0003803404.tif 000032