ライブ蛍光透視画像を用いた冠動脈モデルの非剛体２Ｄ／３Ｄレジストレーション
【課題】心臓インターベンション中の２Ｄ蛍光透視画像を用いたデジタル３Ｄ冠動脈モデルの非剛体レジストレーション方法を提供する。
【解決手段】冠動脈ツリーのデジタル化された３Ｄ中心線の表示を提供するステップ９０と、前記３Ｄ中心線を、少なくとも２つの２Ｄ蛍光透視画像に、グローバル・アライメントするステップ９１と、エネルギー汎関数を最小化することによって、前記３Ｄ中心線を、少なくとも２つの２Ｄ蛍光透視画像に、非剛体レジストレーションするステップ９３と、を含み、前記エネルギー汎関数は、再構築された中心線点とレジストレーションされた中心線点との間の二乗差の総和と、二乗３Ｄレジストレーション・ベクトルの総和と、二乗微分係数３Ｄレジストレーション・ベクトルの総和と、心筋分岐エネルギーと、を含む。
【解決手段】冠動脈ツリーのデジタル化された３Ｄ中心線の表示を提供するステップ９０と、前記３Ｄ中心線を、少なくとも２つの２Ｄ蛍光透視画像に、グローバル・アライメントするステップ９１と、エネルギー汎関数を最小化することによって、前記３Ｄ中心線を、少なくとも２つの２Ｄ蛍光透視画像に、非剛体レジストレーションするステップ９３と、を含み、前記エネルギー汎関数は、再構築された中心線点とレジストレーションされた中心線点との間の二乗差の総和と、二乗３Ｄレジストレーション・ベクトルの総和と、二乗微分係数３Ｄレジストレーション・ベクトルの総和と、心筋分岐エネルギーと、を含む。
【発明の詳細な説明】
【技術分野】
【０００１】
米国出願の相互参照
この出願は、２０１１年９月２８日にハリ・サンダーにより出願された米国仮出願第６１/５４０，１３４号の「心臓インターベンションのガイダンス用のライブ蛍光透視法を用いた冠動脈モデルの非剛体２Ｄ／３Ｄレジストレーション」の優先権を主張するものである。
【０００２】
本発明は、デジタル医用画像の解剖モデルの非剛体２Ｄ／３Ｄレジストレーションの方法に関するものである。
【背景技術】
【０００３】
Ｘ線透視法は、多数の低侵襲的処置、例えば、慢性完全閉塞（ＣＴＯ）の経皮冠動脈インターベンション（ＰＣＩ）のガイダンス用に選択されるモダリティである。これらの処置の間、ガイドワイヤを用いたＣＴＯの横断（crossing）は特に有害である。なぜなら、閉塞が造影剤の伝播を阻止し、蛍光透視法下において、血管の閉塞部を不可視にするからである。さらに、蛍光透視法の射影性質は、３Ｄ構造の解釈を不明確にし、そのことがさらにインターベンションを複雑にする。同様の問題は、低侵襲的処置にも無数に存在している。
【０００４】
現在のプレインターベンショナル計画は、通常、ＣＴ血管造影（ＣＴＡ）あるいは他の３Ｄ画像モダリティの取得を含み、視覚的な不確実性を低減している。血液プール造影剤（contrast）注入を用いたＣＴＡの性能特性は、直接の造影剤注入を用いたインターベンショナルＸ線透視法によって認識可能なものと対照的に、ＣＴＯを引き起こす石灰化を明らかに区別可能にする。これは、計画段階にてＣＴＡを非常に有益なツールとするが、これらのデータとインターベンショナル画像との間の対応を確立することは挑戦すべき問題であると判明している。
【０００５】
この問題に対処するために、処理前の３次元モデルは取得されたボリュームから抽出され、２Ｄ蛍光透視図とアライメントされ、ライブ画像上に重ねられ、それによって、インターベンショナル画像を増加させることができる。しかしながら、このアライメントを達成することは本来困難である。なぜなら、主には、インターモダル対応を見つけることは些細ではないためであり、根本的な最適化の非線形態様のためでもある。さらに、満足な２Ｄ／３Ｄレジストレーションを提供するために、単純な剛体変換では十分でないおそれがある。３Ｄ計画画像が呼吸停止（breathold）の状態で得られるので、自由な呼吸状態で得られる術中の画像と比較して、著しい形状変化が存在する。これは、術前モデルを術中のライブ画像にアライメントする間、非剛体レジストレーション方法を使用することを役立つようにする。しかしながら、これは、困難な設定不備の逆の作業である。それにもかかわらず、２Ｄ／３Ｄレジストレーション方法には、既存の臨床の流れに最小の変更を加えるだけで、インターベンショナルＸ線血管造影法に関する不確定性を大きく低減する可能性がある。より広い展望から、２Ｄ／３Ｄレジストレーション方法は、例えば神経学、整形外科および心臓学の分野における多数の用途を有する。
【０００６】
２Ｄ／３Ｄレジストレーションの最適化面は、画像の離散化、レジストレーションされる構造の複雑さおよび剛体あるいはアフィン変換モデルの使用といったいくつかの理由のために非常に非線形になりうる。その結果、コスト関数の最小化は、困難である。レジストレーションプロセスによって、潜在的に達成可能な最大精度は、含まれる変換モデルの複雑さに直接関連している。多数の２Ｄ／３Ｄレジストレーション方法は、剛体変換モデルあるいはわずかにフレキシブルなアフィン・モデルのみを考慮している。一般的に、これは、骨のような剛体構造や、モダリティの第一のアライメントを与えるには適切であるが、フレキシブルな構造の形状変化を説明するには不十分であると判明している。したがって、非剛体変形モデルは、２Ｄ／３Ｄレジストレーションに対して提案されてきた。
【０００７】
臨床状況において、１つまたは２つの透視面は、インターベンション・ガイダンスのために用いられる。モノプレン（Monoplane）のシナリオが研究されてきたが、報告された誤差は、面外方向において大きい。この正確性の欠如は、レジストレーションプロセスにおける信頼性を著しく制限する。２方向の取得は、これらのインターベンションに関連した不明確性を大きく低減し、今日、日常的に複雑な心臓インターベンションに実行されている。
【発明の概要】
【０００８】
本願明細書に記載されている本発明の実施形態には、概して、２Ｄ／３Ｄ非剛体レジストレーションためのシステムと方法であって、術前の冠動脈の３Ｄ中心線モデルを術中の２以上の蛍光透視画像に非剛体でアライメントし、インターベンショナル画像を視覚的に増加させるシステムと方法が含まれる。本発明の実施形態によるレジストレーションされたモデルは、インターベンショナル血管造影図の上に重ねられ（オーバーレイ可能）、画像導入の慢性完全閉閉塞処理の間、外科的援助を提供し、２Ｄインターベンショナル画像に固有な不確定性を低減する。本発明の一実施形態による方法は、２つの部分、すなわち、（１）第一の固体またはアフィンアライメントを提供するために算出されるグローバル変換モデルと、（２）最終的なレジストレーションを計算するために使用される完全非剛体モデルと、を含む。グローバル・アライメント・パラメータを推定するために、距離マップを用いて、３次元モデルの射影と２Ｄ画像から自動抽出される特徴との間の相違を測定する。第二の部分では、完全非剛体レジストレーション方法を用いて、任意の局所的な形状相違を補償する。本発明の一実施形態による非剛体レジストレーション方法は、変化のフレームワークに基づき、同時照合および再構築プロセスを使用して、非剛体レジストレーションを計算する。最新のハードウェア上の３秒未満の一般的なランタイムでは、本発明の一実施形態によるアルゴリズムは、ＣＴＯ手順の間、手術中に使われるのに十分急速に解析可能である。さらに、本発明の一実施形態による２Ｄ／３Ｄレジストレーション方法を用いて、フレームシーケンス全体にわたって、対象構造の動きをキャプチャすることもできる。
【課題を解決するための手段】
【０００９】
本発明の一実施形態によれば、心臓インターベンション中の２Ｄ蛍光透視画像を用いたデジタル３Ｄ冠動脈モデルの非剛体レジストレーションのための方法が提供され、前記方法は、Qsの３Ｄ制御点からなる一組のSセグメントを具える３Ｄ中心線であって、冠動脈ツリーのデジタル化された３Ｄ中心線の表示を提供するステップと、前記３Ｄ中心線を、少なくとも２つの２Ｄ蛍光透視画像に、グローバル・アライメントするステップと、エネルギー汎関数を最小化することによって、前記３Ｄ中心線を、少なくとも２つの２Ｄ蛍光透視画像に、非剛体レジストレーションするステップと、を含み、前記エネルギー汎関数は、再構築された中心線点とレジストレーションされた中心線点との間の二乗差の総和と、二乗３Ｄレジストレーション・ベクトルの総和と、二乗微分係数３Ｄレジストレーション・ベクトルの総和と、心筋分岐エネルギーと、を含み、前記３Ｄ中心線の非剛体レジストレーションは、前記３Ｄ中心線の座標系にて、対応する中心線点xs,qに適用される一組の３Ｄ並進ベクトルrs,qとして表される。
【００１０】
本発明のさらに別の一実施形態では、前記冠動脈ツリーは前記２Ｄ蛍光透視画像からセグメントされ、前記セグメントされた冠動脈ツリーは前記２Ｄ蛍光透視画像の各々の２値化として表される。
【００１１】
本発明のさらに別の一実施形態では、前記３Ｄ中心線を、前記少なくとも２つの２Ｄ蛍光透視画像に、グローバル・アライメントするステップは、前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線の並進（translation）のみのレジストレーションを実行するステップと、前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線の剛体レジストレーションを実行するステップと、前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線のアフィン・レジストレーションを実行するステップと、
を含み、各ステップでは、
【数１】
の形のエネルギー汎関数は最小化され、ｊは、前記並進のみのレジストレーションと、前記剛体レジストレーションと、前記アフィン・レジストレーションと、を分類し、Cは、前記３Ｄ中心線の座標系と３Ｄ画像形成装置の座標系との間の座標変換であり、Tjは、前記並進のみのレジストレーションと、前記剛体レジストレーションと、前記アフィン・レジストレーションと、に対して、前記３Ｄ画像形成装置の座標系と、前記蛍光透視画像を取得する２方向Ｃアームの座標系と、の間の座標変換であり、xは、均一な座標の３Ｄ点であり、nは、前記蛍光透視画像上の合計を表し、sおよびqは、前記３Ｄ中心線の前記セグメント上の制御点を分類し、Pnは、前記２方向Ｃアームの前記座標系と、第ｎ番目の蛍光透視画像の前記座標系と、の間の座標変換であり、Ψnは、前記３Ｄ位置ｘを前記第ｎ番目の蛍光透視画像に射影する射影オペレータであり、
【数２】
は、前記第ｎ番目の蛍光透視画像の２Ｄ点と、前記第ｎ番目の二値化された画像Bnの前記セグメントされた冠動脈ツリーに最も近い点と、の間の距離である。
【００１２】
本発明のさらに別の一実施形態では、変換行列Tjは、再センタリングされ、前記変換行列Tjの回転パラメータの変化は、前記３Ｄ画像形成装置の前記座標系の原点とは対照的に、前記３Ｄ中心線の前記座標系の原点の周りで発生するように見える運動を誘導する。
【００１３】
本発明のさらに別の一実施形態では、時間的に連続したF対の２Ｄ蛍光透視画像を提供するステップと、前記エネルギー汎関数
【数３】
を最小化することによって、前記３Ｄ中心線を、各対の２Ｄ蛍光透視画像にグローバル・アライメントするステップと、をさらに含み、fは時点の合計であり、Tfは各時点の変換であり、γは所定のパラメータであり、D(Tf・C・xs,q,Tf+1・C・xs,q)はフレーム間距離である。
【００１４】
本発明のさらに別の一実施形態では、前記３Ｄ中心線を少なくとも２つの２Ｄ蛍光透視画像にグローバル・アライメントするステップの後、前記３Ｄ中心線の点を、前記２Ｄ蛍光透視画像の対応する血管中心線点に照合し、前記照合した２Ｄ点から３Ｄ中心線点
【数４】
を再構築するステップをさらに含む。
【００１５】
本発明のさらに別の一実施形態では、前記点照合は、各３Ｄ中心線点xs,pの射影x'を、前記グローバル・アライメントからの各２Ｄ蛍光透視画像上に見つけるステップと、x'から開始して、各２Ｄ蛍光透視画像に対応するセグメンテーション画像Bnを、前記セグメンテーション画像B上の前記３Ｄ中心線の射影に垂直な線に沿って検索し、中心線点xs,pに照合する２つまでの点を見つけるステップと、xs,pに最も近い点を照合点xnとして選択するステップと、前記蛍光透視画像の照合点を用いて、前記心臓インターベンションに用いられている画像処理システムの射影形状を使用している３Ｄ中心線点
【数５】
を再構築するステップと、を含む。
【００１６】
本発明のさらに別の一実施形態では、再構築された中心線点
【数６】
とレジストレーションされた中心線点
【数７】
との二乗差の前記総和は、
【数８】
であり、
【数９】
は指標関数
【数１０】
であり、照合点がいかなる画像平面でも見つからない場合、あるいは、２Ｄ蛍光透視画像上の
【数１１】
の射影x'と前記第ｎ番目の２Ｄ蛍光透視画像の照合点xn"との間の２Ｄ距離が所定の閾値Max2Dより大きい場合、あるいは、中心線点xs,qと対応する再構築された点
【数１２】
との間の３Ｄ距離が所定の閾値Max3Dより大きい場合、H(xs,p)=0であり、その他の場合、H(xs,p)=1である。
【００１７】
本発明のさらに別の一実施形態では、二乗３Ｄレジストレーション・ベクトルの前記総和と、二乗微分係数３Ｄレジストレーション・ベクトルの前記総和と、前記エネルギー汎関数の前記心筋分岐エネルギーは、
【数１３】
であり、rは、レジストレーション・ベクトルであり、xは、前記３Ｄ中心線上の制御点であり、sおよびqの前記総和は、前記３Ｄ中心線の全セグメントの全制御点の合計であり、Kは、前記３Ｄ中心線の接合部から距離Ds離れた各セグメントの点を、前記接合部で接合される他のセグメントの対応する制御点に連結することにより形成される制御点の全対の一組であり、前記接合部では、連結がレジストレーションの間拡大または圧縮される場合、各連結はその変位に比例した力を発生させる。
【００１８】
本発明のさらに別の一実施形態では、前記エネルギー汎関数は、
【数１４】
を反復することによって、最小化され、zは、点rs,p=[0,0,0]T|∀{s,q}から反復z=0で開始する反復カウンタであり、εは、小さい定数であり、Ks,q⊂Kは、制御点(s,q)に作用する一組の制約であり、μ、υおよびλは、所定のパラメータである。
【００１９】
本発明のさらに別の一実施形態では、コンピュータ可読の非一時的プログラム記憶装置が提供され、前記コンピュータにより実行されるとともに、心臓インターベンション中の２Ｄ蛍光透視画像を用いたデジタル３Ｄ冠動脈モデルの非剛体レジストレーションのための方法を実行する命令プログラムを有する。
【図面の簡単な説明】
【００２０】
【図１】本発明の一実施形態に係る、対象とするシステムの形状の概略図である。
【図２】（ａ）〜（ｈ）は、発明の一実施形態に係る、レジストレーションの複数のステージによる入力データの変換を示す図である。
【図３】発明の一実施形態に係る、非剛体レジストレーションの進行を示す図である。
【図４】発明の一実施形態に係る、中心線制御点照合アルゴリズムのフローチャートである。
【図５】（ａ）（ｂ）は、本発明の一実施形態に係る、点照合および分岐点における分岐制約の形成を示す。
【図６】本発明の一実施形態に係る、反復回数の関数として、中心線の全移動を示すグラフである。
【図７】本発明の一実施形態に係る非剛体レジストレーション方法において、パラメータλの挙動を示す図である。
【図８】本発明の一実施形態に係る非剛体レジストレーション法の後、平均２Ｄおよび３Ｄ誤差を示す図である。
【図９】本発明の一実施形態に係る、非剛体２Ｄ／３Ｄレジストレーション法のフローチャートである。
【図１０】（ａ）〜（ｄ）は、本発明の一実施形態に係るＤＤＲの形成を示す図である。
【図１１】本発明の一実施形態に係る初期位置の摂動に関して、オプティマイザの性能を例示する。
【図１２】本発明の実施形態に係るノイズ・レベルに対する入力ＤＲＲノイズの標準偏差に関して、オプティマイザの性能を例示する。
【図１３】本発明の実施形態に係るシミュレーションデータを有する若干のサンプル非剛体レジストレーションを示す。
【図１４】本発明の実施形態に係るシミュレーションの非剛体変形レベルに関する３Ｄ残差エラーのグラフである。
【図１５】本発明の実施形態に係る、剛体、アフィンおよび非剛体レジストレーション後の入力画像および結果を表す。
【図１６】本発明の実施形態に係る、５人の患者のための剛体、アフィンおよび非剛体レジストレーション後の２Ｄマニュアル・セグメンテーションに関して２つの透視面上の中心線の平均残余の２Ｄ射影エラーを示している表である。
【図１７】本発明の実施形態に係る、一連のフレーム上のレジストレーションを示す。
【図１８】本発明の実施形態に係る、１つの心周期の複数のフレームによるＲＣＡのトラッキングを表す。
【図１９】本発明の一実施形態に係る、ライブ蛍光透視画像を用いた冠動脈モデルの非剛体２Ｄ／３Ｄレジストレーションのための方法を実施するためのコンピュータシステムの一例のブロック図である。
【発明を実施するための形態】
【００２１】
本発明の例示的実施形態は、概して、ライブ蛍光透視画像を用いた冠動脈モデルの非剛体２Ｄ／３Ｄレジストレーションのためのシステムおよび方法を含む。したがって、本発明は各種変更および代替実施形態の影響を受けるが、特定の実施形態を、単に一例として図面および明細書に詳細に記載する。しかしながら、本発明は、記載されたこの特定の形態に限定することを意図するものではなく、本発明は、本発明の精神および範囲内にある全ての変更、等価および代替を含むものである。
【００２２】
本明細書に記載の用語「画像」は、個別の画像要素（例えば、２次元画像のピクセルおよび３次元画像のボクセル）から構成された多次元データを意味する。画像は、例えば、コンピュータ断層撮影法、磁気共鳴映像法、超音波、あるいは、従来当業者に知られている他の任意の医用画像システムによって収集された被写体（subject）の医用画像とすることができる。画像はまた、医学とは関係ない分野、例えば、リモートセンシングシステム、電子顕微鏡法等から提供される場合もある。画像は、R3からRまたはR7の関数とみなすことができるが、本発明の方法はこれらの画像に限定されず、任意の次元、たとえば、２次元画像あるいは３次元体積の画像に適用可能である。２次元あるいは３次元画像のために、画像の定義域は一般的に２次元あるいは３次元の長方形配列であり、各ピクセルあるいはボクセルは、２本あるいは３本の相互に直交する軸に関して方向付け可能である。本明細書に記載の用語「デジタル」および「デジタル化」は、デジタル取得システムを介して、あるいはアナログ画像からの変換によってデジタル化されたフォーマットにおける画像あるいは体積を適切に意味する。
【００２３】
本発明の一実施形態によると、術前のＣＴＡからセグメントされた、冠動脈ツリーの３Ｄ中心線表現は、２以上の同時の蛍光透視画像に非剛体レジストレーション可能である。初期設定の３Ｄ位置から開始し、装置のキャリブレーションから計算し、中心線の３Ｄアライメントは、並進のみの運動と剛体変換とアフィン変換とを用いて徐々に調節される。最後に、非剛体変換が計算され、最終的なレジストレーションが提供される。
【００２４】
本発明の一実施形態によると、冠動脈ツリーは、制御点xs,p∈R、s∈{1,...,S}、q∈{1,...,QS}から構成された一組のSセグメントQsによって表され、方向付けのない（undirected）非巡回グラフを形成する。説明を明確にし、一般性を失わないために、以下の説明は、１対の蛍光透視画像In、n∈{1,2}が与えられる場合に制限されている。しかしながら、この制限は例示であり限定的ではなく、本発明の実施形態に係る方法は、複数の蛍光透視画像に拡張可能であることを理解されたい。さらに、説明を明確にし、一般性を失わないために、両者が同一のサイズW×H∈N2であると仮定する。血管は、既知の方法を用いて入力蛍光透視画像から自動的にセグメント化され、入力画像In(i,j)の２値化Bn(i,j)として表される。位置(i,j)∈{0,...,W-1}×{0,...,H-1}における入力画像ピクセルが対象の構造に対応する場合、Bn(i,j)=1であり、その他の場合、Bn(i,j)=0である。
【００２５】
図１には、対象とするシステムの形状が概略的に示され、５つの座標系（ＣＳ）を用いて記載されている。図１中のラベルつきのフレームは座標系を表し、破線および太字の符号は変換行列を表す。冠動脈ツリーを表す３Ｄ中心線はＵに関して記載され、座標系（ＣＳ）Ｃは３Ｄ画像装置の基準点に中心がある。２Ｄ血管造影法に関して、Ｘは２方向Ｃアームの基準ＣＳであり、ＣＳのP1およびP2は、２Ｄ結像面に中心を置いている。ＣとＵとの間の変換およびＸと{P1,P2}との間の変換は、装置のキャリブレーションから知られており、それぞれ、剛体変換行列C、P1およびP2で符号化される。２つのＸ線平面の射影形状もまた利用可能に仮定され、射影オペレータΨ1およびΨ2で符号化される。射影オペレータΨ1およびΨ2は、同次座標x=[x,y,z,1]Tの３Ｄ点を、対応する透視面上の２Ｄ点x’=[x,y]にマップする。２つの画像診断座標系の間の変換は、Tによって表される。
【００２６】
完全性のために、例示であり限定的ではないが、自動セグメンテーション方法を以下に簡単に説明する。（１）入力蛍光透視画像は、ヘシアンベースのサトー血管フィルタ（Hessian-based Sato vesselness filter）を用いて処理され、（２）血管測定の閾値化をまばらにサンプリングすることによって、シードポイントが生成され、（３）これらのポイントから開始し、各ピクセル位置にて計算されたヘシアン行列の第１の固有ベクトルに沿って統合されることによって、ファイバが生成され、（４）ファイババンドルは反復して間引かれ（thinned）、血管の中心線を引き出す。
【００２７】
並進、剛体およびアフィン変換
一般的な臨床状況では、３Ｄ冠動脈中心線をＣＴＡ取得から２方向透視における同一構造の画像にアライメントする変換は、大まかにしかできない。なぜなら、患者の位置は高い精度で制御することができないためである。本発明の一実施形態に係る変換T、すなわち、２つの画像診断座標系間のマッピングは、装置のキャリブレーションによって推定可能であり、対象の構造の形状を考慮することによって調整可能である。初期条件から開始し、改良されたアライメント変換（alignment transformations）は、最小化プロセスを用いて計算され、グローバル変換モデルあるいは剛体変換モデルのパラメータが推定される。この点から、アフィン変換モデルもまた用いられ、３Ｄ中心線を変形し、ＣＴＡと２方向Ｘ線との間の形状の不一致を補償する。
【００２８】
本発明の一実施形態によると、３Ｄ中心線の射影と、各透視面からセグメント化された血管中心線との間の総距離は最小化される。χ(x)=T・C・xは、Ｕに関する点x=[x,y,x,1]をX座標系にマップするアフィン変換演算子とする。また、Пn(x)=Ψn・Pn・xは、Xからの点を透視面Pnにマップする射影である。本発明の実施形態は、エネルギーを、２つの蛍光透視画像で最小化され加算されるべき量として、以下のように定義する。
【数１５】
ここで、
【数１６】
は２Ｄ点とBn(i,j)=1である最接近点との間の距離である。各画像位置に対する距離
【数１７】
は、各Bn(i,j)セグメンテーションの距離変換を形成することによって、事前に計算される。２つの蛍光透視画像を使用することは例示であり限定ではなく、２つより多い蛍光透視画像が用いられる場合、加算は容易に拡張されることを理解されたい。
【００２９】
図２（ａ）〜（ｈ）は、入力データの変換を示している。図２（ａ）はＣＴ体積からの湾曲した（curved）スライスを示し、図２（ｂ）は３Ｄ中心線を示し、図２（ｃ）（ｄ）は、異なる視点からの透視血管造影図を示し、図２（ｅ）（ｆ）は蛍光透視画像のセグメンテーションを示し、図２（ｇ）（ｈ）はセグメンテーションの距離変換を示す。
【００３０】
Tの１２の値を直接最適化することは、効果的でなく非効率なため、本発明の実施形態では、代わりに、パラメトリック手法を用いる。並進のみの変換TT、剛体変換TRおよびアフィン変換TAは、それぞれ、パラメータセットT∈R3、R∈R6およびA∈R12を用いて表される。TSのパラメトリック表現から行列表現へのマッピングは、以下のように定義される。
【数１８】
【数１９】
【数２０】
【００３１】
ここで、Sijは、Sの添え字{i, ... ,j}を伴う要素を含む部分集合である。
【数２１】
【数２２】
【数２３】
ここで、cx=cos(θx)、sx=sin(θx)等である。本発明の実施形態では、変換行列TT、TRおよびTAは、再センタリングされる。回転パラメータのいかなる変化も、Ｃの原点と対照的に、Ｕ座標系の原点周辺で発生するように見える運動を誘導するというような方法で、この動作は、式（１）のエネルギー関数に影響しないが、通常オプティマイザに有益である。なぜならこれは、直観的に、ＣＴボリュームの重心で発生するように見える回転動作またはスケーリング動作に関連して、この種の再センタリングが、パラメータ空間において、距離を減らすためである。最終的なTの定義は、式（４）のとおりである。
【数２４】
ここで、ＣはＵからＣへの変換であり、j∈{T,R,A}であり、選択された変換モデルに依存する。式（１）は、３つの変換モデルの各々に対して、式（１ａ）の形をとなる。
【数２５】
任意の非派生的な多目的最適化アルゴリズムを用いて、式（１）を最小化できる。
【００３２】
初期設定の患者の位置から開始して、本発明の実施形態は、（１）並進のみの変換、（２）剛体変換および（３）アフィン変換を考慮することによって、アライメントを徐々に調整する。
【００３３】
マルチ・フレーム・アラインメント
本発明の一実施形態に係るグローバル・アライメント方法は、マルチ・フレーム・シナリオをカバーするために拡大されることができる。この場合、実施形態は、時間的連続性のため、推定されたアライメント変換が１フレームからその他にほとんど何も変えないと仮定する。本発明の一実施形態に係るグローバルなエネルギー関数は、一連の２方向フレームのために以下のように、再定義可能である。
【数２６】
あるいは
【数２７】
である。ここで、γは自由パラメータであり、χfは、χと同様に定められるが、時点ごとの異なる変換行列Tfを有する。変換行列Tfは、上述したように、並進のみ、リジッドあるいはアフィンとすることができる。上述したように、２つの蛍光透視画像上の合計nは例示であり限定的ではなく、２つ以上の蛍光透視画像が使われている場合、合計は実施形態において、容易に拡張されうる。
【００３４】
フレーム間距離D(χf,χf+1)は、リジッド変換あるいはアフィン変換とすることもできる。リジッド変換行列およびアフィン変換行列の空間が線形でないので、本発明の一実施形態に係るフレーム間距離D(χf,χf+1)の定義は精巧である。本発明の実施形態により仮定されているように、剛体変換空間は方向における小さい変化のために局所的に線形化可能である。このフレームを用いて、フレームfからフレームf+1への剛体変換は、並進ベクトル
【数２８】
および回転ベクトル
【数２９】
によって表され、
【数３０】
となるように、単位ベクトル
【数３１】
および回転角度θf,f+1のプロダクトとして定義される。並進ベクトルは、剛体変換行列
【数３２】
の並進部分に対応し、行列
【数３３】
と回転ベクトル
【数３４】
との間のマッピングは、
【数３５】
および
【数３６】
によって与えられる。
【００３５】
この表現を使用して、実施形態は、以下の通りにフレーム間距離を形式的に定めることができる。
【数３７】
【数３８】
ここで、αは、回転部分および並進部分のそれぞれの貢献のバランスをとる実数である。本発明の実施形態は、その内容全体が本願明細書に参照によって取り込まれる、ボワヴェール（Boisvert）等による「Geometric Variability Of The Scoliotic Spine Using Statistics On Articulated Shape Models」（ＩＥＥＥ Ｔｒａｎｓ. Ｍｅｄ. Ｉｍａｇ.、第２７巻、第４号、５５７〜５６８頁、２００８年４月）に従って、αを０．０５に設定する。
【００３６】
アフィン変換行列の間の線形距離を定めることは、関連が深い複雑な計算を必要とする。本発明の実施形態では、剛体変換が変換の大部分を捕えると仮定されるので、本発明の実施形態は、変換の剛体部分を正規化し、残りのアフィン・パラメータを無拘束のままにする。形式的に、アフィン変換の場合、正規化定義の変換TfはパラメータA16を使用して構成されるが、１２のパラメータのフルセットはエネルギー定義の変換χfを構成するために用いられる。
【００３７】
このエネルギー関数は、フレームシーケンスの上のアライメント変換をトラッキングするのに役立ち、または、γ→∞の場合、例えば、３つの隣接するフレームを使用して平均的変換を計算するのに役立つ。シングル・フレーム・エネルギーと同様に、変換χfは、変換モデルの複雑さを順次増加させることによって、次第に精製される。
【００３８】
非剛体レジストレーション
グローバル・アフィン変換モデルが、呼吸によって、および、心臓の鼓動によって、生じるインターモダル形状矛盾を完全に補償できるというわけではない。本発明の実施形態は、３Ｄ冠状動脈ツリー中心線と、２つのキャリブレーションされた蛍光透視画像と、の間の視覚的な相関関係を改善できる非剛体レジストレーション方法を提供する。本発明の一実施形態に係る方法は、ユーザが相関関係を識別する必要がないという点で自動といえる。
【００３９】
本発明の一実施形態に係る非剛体変換は、３Ｄ並進ベクトルrs,qのセットrとして表され、３Ｄ並進ベクトルrs,qは、座標（ＣＳ）のＵにおける対応する中心線点xs,pに与えられる。次に、レジストレーション点
【数３９】
は
【数４０】
を用いて計算される。溶液の品質を測定するために規定された、本発明の一実施形態に係るエネルギーは、以下（８）式の形を有する。
【数４１】
本発明の一実施形態に係る最小化プロセスは反復され、位置rs,p=[0,0,0]T|∀{s,q}から開始し、変換rs,qは、勾配降下法に従うことによって徐々に調整される。これは、レジストレーションされた中心線の滑らかな動きを結果として生ずる。図３は、非剛体レジストレーションの進行を示し、図中左から右に、０回、３５０回および１２００回の反復後のサンプル入力画像および中心線位置を示している。
【００４０】
画像エネルギー
画像エネルギーという用語は、再構築された中心線の汎関数である。本発明の一実施形態によると、照合および再構築プロセスを用いて、３Ｄ中心線モデルにおいて全ての点xsqに対する再構築点
【数４２】
を計算することによって、各蛍光透視画像における中心線の射影を構成することができる。図４のフローチャートを参照して、照合は各透視面上で別個に以下のように行われる。本発明の一実施形態に係る照合方法は、ステップ４１において、３Ｄ中心線の制御点を蛍光透視画像nに射影することによって開始する。蛍光透視画像n上の３Ｄ中心線点
【数４３】
の射影x′から開始して、ステップ４２において、蛍光透視画像nに対応するセグメンテーション画像Bn上の検索が、中心線の射影に垂直な方向に行われる。この方向は、各位置において、偏導関数
【数４４】
に基づいて計算される。その結果、ステップ４３において、多くても２つの照合点xa"、xb"は各透視面n∈{1,2}上で識別される。ステップ４１〜４３は、他の蛍光透視画像のためにステップ４４から繰り返される。ステップ４５において、エピポーラ拘束
【数４５】
を最小化するペア(x1,i",x2,i"),i∈a,bが保たれる。ここで、d(x1,i",em,i")は点xn,i"と、xm,i"によって導入されるエピポーラ線との間のユークリッド距離であり、{n,m}∈{1,2}かつn≠mである。ステップ４６において、式（９）を最小化するペアは、システムの射影形状を使用して、再構築３Ｄ点
【数４６】
を計算するために用いられる。ただ１つの点xn,i"が画像面上に見られる場合、それは初期設定により保たれる。図５（ａ）は、中心線x'からの点と２Ｄ構造からの点との照合を示している。点xa"、xb"は、有力候補として識別され、xb"はx'に近いため保たれる。次に、ステップ４１〜４６は、３Ｄ中心線の全制御点のためにステップ４７から繰り返される。
【００４１】
本発明の実施形態では、中心線の射影に局所的に垂直な方向において照合点の検索が行われる。構造のグローバル・アライメントが正しいということは先験的に知られているので、垂直検索方向は、最短点手法と比較して、多数の非現実的な照合を回避する。このことは、２Ｄ自動セグメンテーション（automatic segmentation）が完璧ではない場合に有効である。例えば、３Ｄ射影が、自動的にセグメント化された２Ｄ曲線より長い曲線になった場合、最短点をとることは、チップからの曲線の著しい短縮である。一方、３Ｄ中心線射影のチップにおいて照合は見出されないため、垂直方向での検索は、この効果を低減する。
【００４２】
再構築の品質は保証されず、本発明の実施形態は指標関数
【数４７】
を用いて異常値を棄却する。（１）あらゆる画像面において照合点が見つからない場合、（２）x'とxn"との間の２Ｄ距離がMax2Dより長い場合、（３）xs,qと
【数４８】
との間の距離がMax3Dより長い場合、H(xs,p)=0である。ここで、Max2DおよびMax3Dは、所定の閾値である。実施形態によると、一例として、限定的ではないが、Max2Dは５０ピクセルであり、Max3Dは６ｍｍである。
【００４３】
実施形態によると、式（８）における画像エネルギーは再構築点
【数４９】
とレジストレーションされた中心線点との間の二乗誤差の合計として、全ての有効な点に関して以下のように定義される。
【数５０】
全てのレジストレーションされた中心線点が再構築点と同一の位置にあるとき、このエネルギーは最小となる。
【００４４】
内部エネルギー
本発明の一実施形態によれば、正規化（regularizer）を用いて、視覚的にレジストレーションされた中心線を視覚的にコヒーレントにかつ幾何学的に妥当に保つことができる。本発明の実施形態に従う以下の３つの内部エネルギーの項は、このために考慮可能である。
【数５１】
ここで、μ、νおよびλがエネルギー・バランス・フリー・パラメータである。第１の項EDisp(r)は、非剛体変換による変位が小さいとき最小である。この正規化は、一対の２Ｄ血管造影によって、わずかにしか拘束されない３Ｄ面外動作を最小化する。第２の項ESmooth(r)は、各血管セグメント上の平滑性を確実にすることができる。ここで、
【数５２】
および
【数５３】
は、それぞれ、セグメントパラメータに対する位置(s,q)のための並進ベクトルｒの第１および第２の微分係数を表す。セグメント接合部で、ノイマン境界条件
【数５４】
が仮定される。
【００４５】
はじめの２つの内部エネルギー項については、中心線モデルはフレキシブルであり、接合ノードにおいて以外は、血管の動作は互いに独立している。実際は、血管動作は、血管自身の固有の剛性によってだけでなく、血管が心筋に取り付けられているという理由によっても機械的に拘束される。この点において、第３の項、EMyocard(r)は、心筋制約の最小モデルとしての機能し、血管分岐で剛性の特定の程度を確保し、小さいセグメントがより大きいセグメント上へ倒れこむのを防止できる。本発明の実施形態は、図５（ｂ）に示すように、３つのセグメントの接合部の周りの人工的なリンクを用いて、この制約をモデル化し、図５（ｂ）は、接合部における３つの心筋制約（k1、k2、k3）の形成を示す。Dsを距離パラメータとし、Nに設定される値と接合部に存在する血管の平均直径とを掛け算する。Nの例示的かつ非限定的な値は３である。接合部から距離Ds程離れた各セグメント上のノードは、機械ばねに類似した方法で、他のセグメントの対応するノードに連結され、各リンクは、レジストレーションプロセスの間、拡大または圧縮される場合、その変位比率と比例する力を発生させている。式（１１）において。Kは、ともに結合された全対のノード{i,j}の組として定義される。
【００４６】
エネルギー最小化
本発明の一実施形態による完全なエネルギー関数は、式（８）、式（１０）および式（１１）を使用して定められ、そのオイラー・ラグランジェ方程式を計算して、勾配降下方法に従うことによって、本発明の実施形態において最小化できる。本発明の一実施形態によれば、反復rs,p=[0,0,0]T|∀{s,q}かつ反復z=0での点から開始して、結果として生じる離散更新方程式は、式（１２）のとおりである。
【数５５】
ここで、数ステップεは、小さい定数であり、Ks,q⊂Kは、ノード(s,q)に作用する制約の中のエンプティセットである。本発明の実施形態では、式（１２）は、全体の変位が収束するまで、反復的して評価されうる。より形式的には、全体の変位は、
【数５６】
を意味し、停止基準は、ｚ＞βに関してΔz-Δz-β<ζΔβとして定められる。図６は、本発明の一実施形態に係る、非剛体レジストレーションの間、反復回数の関数として、中心線の全体の変位を示している。曲線の形状から推論されうるように、最終結果はβおよびζの特定の値にあまり影響されない。例示的かつ非限定的な値は、β＝100およびζ＝0.001である。
【００４７】
パラメータ選択
本発明の一実施形態によれば、３つのパラメータμ、υおよびλは、式（１２）の挙動を支配し、特定のタスクに対して選択可能である。特定のパラメータセット（μ＝0.10、υ＝10.0、λ＝1.0）から開始して、以下に詳細に記載されているように、各値を個別に変化させることによって、平均２Ｄ射影エラーおよび３Ｄエラーを計算することによって、方法の挙動についての洞察が得られる。この方法は、後述するように、アフィンアライメントの後、臨床データセットの患者４、および、合成変形を有するデータセットに適用された。前者はＲＣＡを表し、後者はＬＣＡを表す。図８は、非剛体レジストレーション後の平均２Ｄおよび３Ｄエラーを示す。図８において、上の行はパラメータμのデータを示し、中間の行はパラメータυのデータを示し、下の行はパラメータγのデータを示す。各グラフに対して、１つのパラメータが変化し、他の２つのパラメータは一定に保たれた。デフォルト・パラメータ値は、μ＝0.1、υ＝10.0、λ＝1.09である。ｘ軸のスケールは、対数関数である。
【００４８】
パラメータμは全体の変位（displacement）を拘束することができ、２Ｄエラー曲線は計算されたエラーがその値とともに増加することを示している。しかしながら、３Ｄエラー曲線は、μが０．０５から０．１０の範囲の値を有するとき、最小値を示す。加えて、収束率は、μとともに増加する。例えば、レジストレーションを用いて、μ＝{0.001,0.01,0.1}をそれぞれ有する４８８０、３３２０および、１１００の反復に収束した図８（左上）のグラフを生成する。
【００４９】
パラメータυは、モデルの剛性を制御できる。この場合、２Ｄ射影エラーは、０．５０から２．００の範囲の値を使用することを提案する。しかしながら、３Ｄエラー曲線は、より高い値がより正確なレジストレーションになりうることを示す。本発明の実施形態は、例示的かつ非限定的な値として、２Ｄエラー曲線と３Ｄエラー曲線との間の良好なトレードオフを示すυ＝10.0を使用する。
【００５０】
パラメータλは、心筋制約のバランスをとり、グローバル形状の一般的な整合性を保存できる。臨床ＲＣＡデータセットについては、高い値が運動を過度に拘束し、エラーの増加を引き起こすことが判明した。しかしながら、シミュレーションされたＬＣＡデータセットについては、状況はほぼ逆に現れ、１０．０までのλ値は、明らかに有益である。この状況は、ＬＣＡおよびＲＣＡのデータセットに影響を及ぼす変形の性質における違いによって、おそらく説明可能である。例示的かつ非限定的な値λ＝0.5が選択された。なぜなら、それがエラー測定のトレードオフを表し、かつ、図７にて示すように、多くのケースの視覚的外観を改善するためである。図７は、非剛体レジストレーション方法のパラメータλの挙動を、左側のλ＝0および右側のλ＝0.3に関して示している。ＬＣＡおよびＲＣＡデータセットに対して、複数の値のλを用いることができる。本発明の実施形態では、特に言及しない限り、以下の全ての実験において、例示であり限定的ではないパラメータμ=0.5、υ=10.0およびλ=0.5が用いられる。
【００５１】
図９には、本発明の一実施形態に係る非剛体２Ｄ／３Ｄレジストレーション方法のフローチャートを示す。方法は、ステップ９０において、術前の画像体積、例えばＣＴＡ体積からセグメント化された冠動脈ツリーの３Ｄ中心線表現を提供することによって開始する。方法は、ステップ９１において、２つ以上の分割された、冠動脈血管の２Ｄ蛍光透視画像に対する冠動脈ツリーの３Ｄ中心線表現の最初のグローバル・アライメントを実行することによって続行する。実施形態では、上述したように、このグローバル・アライメントは、連続的な並進のみの変換、剛体変換およびアフィン変換によって実行可能である。ステップ９２において、図４を参照して説明したように、３Ｄ中心線の点は、蛍光透視画像における対応する血管中心線点に照合され、照合された２Ｄ点から３Ｄ中心線を再構築する。ステップ９３において、非剛体レジストレーションを実行する。次に、ステップ９３において、エネルギー汎関数を最小化することによって非剛体レジストレーションを実行する。ここで、エネルギー汎関数は、再構築点
【数５７】
とレジストレーションされた中心線点との間の二乗差の総和と、二乗３Ｄレジストレーション・ベクトルの総和と、二乗微分係数３Ｄレジストレーション・ベクトルの総和と、心筋分岐エネルギーと、を含む。
【００５２】
実験結果
本発明の実施形態に係るアルゴリズムは、インターベンション中にレジストレーションされたモデルのライブ蛍光透視画像に対するオーバレイを示すプロトタイピング環境においてＣ＋＋で実施された。Ｃ＋＋での実施は例示であり限定的ではなく、本発明の他の実施形態では、他のコンピュータ言語での実施が可能であることを理解されたい。全てのグローバル・オプティマイザはシングルスレッドである。本発明の一実施形態に係る非剛体レジストレーション法は、OpenMPを使用して部分的にマルチスレッドされている。与えられるランタイムはi7 Q820 Quad Core Intel CPU用である。示されたランタイムは２Ｄ画像セグメンテーションを構成しない。このオペレーションは、画像ごとに平均して０．２００秒のランタイムを有し、標準偏差は０．０８８秒である。シミュレーションと臨床データの両方に関する実験が行われた。
【００５３】
シミュレーション
本発明の実施形態によるアライメントおよびレジストレーションの方法は、シミュレーションのセットにおいてテストされ、いくらか理想的なシナリオにおける本発明の実施形態に従う各種アルゴリズムの性能特性を与え、臨床データを使用する場合には不可能な３Ｄレジストレーションエラーを評価した。中心線は、以下の通り生成され、５１２×５１２ピクセルのサイズおよび０．３４５ピクセル／ｍｍの解像度を有する、一対のデジタル的に再構築されたＸ線画像（ＤＲＲ）にレジストレーションされた。（１）ＣＴボリュームにおいて、セグメントされた冠動脈は、シミュレーションされた透視面に投影され、（２）結果として生じる射影の逆数は、参照背景画像から減算された。図１０（ａ）〜（ｄ）は、ＤＲＲの作成を例示する。ＣＴボリュームにおいて、セグメントされた冠状動脈血栓は、図１０（ａ）のシミュレーションされたＸ線撮影平面に投影される。次に、この画像は、図１０（ｂ）に示すように実際の造影剤のない（contrast-free）血管造影法から減算され、図１０（ｃ）の最終的なＤＲＲを生成する。一実験において、図１０（ｄ）に示すように、σ＝50を有するガウス雑音が加えられる。シミュレーションは、後述する３つのシナリオにおいて使われた。すべての変換がシミュレーション中に公知であるので、本願明細書において提示されるすべての３Ｄエラーは、正確なポイント・ツー・ポイント・エラーである。視覚効果を生じないときでも、エラーのこのスタイルがセグメントの圧縮または拡大を不利にすることに留意する必要がある。それゆえ、それは、大部分のＴＲＥ形成より厳格である。２Ｄエラーは、上述したように、参照セグメンテーションとして実際の中心線の射影を用いて計算される。
【００５４】
（１）初期の解決法に対する依存
中心線は一対のＤＲＲによってレジストレーションされ、異なる初期位置から開始する。これは、初期位置の摂動と、同様に、C、P1、P2またはTの剛体変換行列のミスキャリブレーションに対するグローバル・アライメント方法の感度を定量化する機能を果たす。選択された初期位置は、以下の通りである。T0を装置のキャリブレーションによって、与えられる変換であると仮定する。次に、１２の初期ポイントは３本の主軸に沿って±αｍｍだけ変位したT0と、は３本の主軸の周りに±β度だけ回転したT0である。実施形態は、α∈10×{0,1,...,6}mmおよびβ∈3.75×{0,1,...,6}度を使用する。平方二乗平均エラーは、摂動（並進または回転）の同じレベルを有する全結果に対して計算された。図１１は、初期位置の摂動に関して、オプティマイザの性能を例示するものであり、左列のグラフは剛体アライメントの結果を示し、右列のグラフはアフィンアライメントの結果を示し、上行は方向における変化（変位）の結果を示し、下行は位置における変化の結果を示す。下行において、失った点は、平方二乗平均エラーが３ｍｍ超であることを表す。図１１の結果は、回転摂動が１０度未満のとき、提示される本発明の実施形態に係るすべてのアルゴリズムが満足であることを示す。
【００５５】
（２）画像ノイズに対するロバスト性
分散を有する複数の量のガウス雑音は、蛍光透視画像に加えられ、本発明の実施形態に従う方法の性能上の画質品質の影響を評価した。各オプティマイザは、各ノイズ・レベルに対する１３の初期位置から開始した（上述したように、α＝5mmおよびβ＝7.5度がプラスされたT0）。図１２は、本発明の実施形態に従って、ノイズ・レベルに対する入力ＤＲＲノイズの標準偏差に関して、オプティマイザの性能を例示し、左列のグラフはリジッドアライメントの結果を表し、右列のグラフはアフィンアライメントの結果を表す。σ∈[0,60]を有するガウス白色雑音モデルが用いられた。この試験の結果は、本発明の実施形態による方法が蛍光透視のガウス雑音の存在に影響されないことを示す傾向がある。臨床設定ではほとんど観察されないσ＝60のレベルでさえ、２Ｄおよび３Ｄエラー・レベルは、追加されたノイズの量にわずかな相関しか示さない。これは、選択された２Ｄセグメンテーション・アルゴリズムの性能が良好であるという徴候である。
【００５６】
（３）非剛体変形
これらのシミュレーションにおいて、中心線は、ＣＴで走査した領域をカバーしている３×３×３ノード薄板スプライン（ＴＰＳ）変形モデルを使用して変形した。各ＴＰＳノードは、変形パラメータξに応じた係数によって、中心ノードへ移動する。９つの左のノードは、係数ξによってシフトされ、９つの中心のノードは、0.5によってシフトされ、９つの右のノードは、係数2ξによってシフトされた。図１３は、シミュレーションされたデータを有する若干のサンプルの非剛体レジストレーションを示す。上行は初期位置を示し、下行は最終位置を示す。本発明の一実施形態による非剛体レジストレーション・アルゴリズムは、変形レベルξ∈{0.020,0.025,...,0.060}に対して試験され、中心線の平均３Ｄ変位は、［1.209〜3.638］ｍｍであり、中心線の最大３Ｄ変位は、［2.472〜7.416］ｍｍであった。非変形中心線曲線を、初期の解決法として用い、非剛体レジストレーション・アルゴリズムを適用する前に、グローバル・アライメントは、実行されなかった。図１４は、定量的評価を示しており、シミュレーションされた非剛体変形レベルに関する３Ｄ残差のグラフであり、ベースラインは初期の解決法のエラーに対応する。このアルゴリズムがレベルξ＝0.045までは良好な結果を示すことが判明した。この点では、平均的３Ｄエラーは２．７２１ｍｍから１．１９８ｍｍに減少し、このことは、図１３に示すように質的に重要である。以下、本願明細書に示されるように、性能のこのレベルは臨床シナリオにおいて適切である。
【００５７】
臨床データ
以下のセクションに示される実験に、５つのデータセットが使われた。各データセットは、拡張末期（データセット１、３）または収縮末期（データセット２、４、５）にて得られる１つのＣＴＡ走査と、１つの２方向Ｘ線透視法記録と、を含む。２つのモダリティはＥＣＧゲーティングを用いて時間的にアライメントされ、冠動脈は専門家によって、ＣＴＡにおいて、半自動式にセグメントされた。セグメントされた中心線の平均３Ｄインターポイント距離は、［1.36〜1.60］ｍｍの範囲にあった。血管造影図は、３つの場合の左冠状動脈（ＬＣＡ）と、その他２の場合の右冠動脈（ＲＣＡ）と、を撮像する。画像サイズは全ての場合において、５１２×５１２であり、画像解像度はデータセット１〜５に対してそれぞれ｛0.345,0.279,0.279,0.216,0.279｝ｍｍ／ピクセルである。獲得率は、１５ｆｐｓである。標準Ｃアームキャリブレーション行列が使われた。必須ではないが、キャリブレーション行列は、データセットごとに異なるにもかかわらず、獲得の間一定であった。ＣＴＡおよび蛍光透視試験は、冠動脈疾患の治療のための患者に定められた。
【００５８】
本発明の一実施形態によるアライメントおよびレジストレーション・アルゴリズムの性能の評価は、平均２Ｄ射影エラー
【数５８】
を使用して定量化された。ここで、Mは３Ｄ中心線モデルの点の総数であり、
【数５９】
は参照２ＤセグメンテーションSnの距離変換であり、対応する蛍光透視画像上で手動追跡によって得られた。
【００５９】
グローバル・アライメント：オプティマイザの性能評価
本発明の実施形態は、グローバル・アライメントのために、７つのローカル（Ｌ）アルゴリズムおよび２つのグローバル（Ｇ）アルゴリズムを検査した。以下、これらのアルゴリズムを簡潔に説明する。
【００６０】
ベスト・ネイバー（Best Neighbor）（Ｌ）
各反復で、パラメータ空間のすべての主要な方向におけるステップが考慮される。最善の移動が適用され、新規の反復が開始する。移動が解決法を改善しないとき、ステップは半分にされる。
【００６１】
ネルダー・ミード（別名滑降シンプレックス）（Ｌ）
一般のシンプレックス（またはポリエイラクブカ）のｍ＋１頂点で関数値を評価することによって、ｍ次元関数を最小化する古典的な数の最適化方法である。各反復で、最悪の値を有する頂点は、拡大動作または縮小動作の一方が続く反転動作を使用して、別のものに置換される。ネルダー等による「A Simplex Method For Function Minimization」、Ｃｏｍｐｕｔ．Ｊ．第７巻、第４号、３０８〜３１３頁、１９６５年参照。
【００６２】
Ｓｂｐｌｘ３（Ｌ）
この方法は、低次元の部分空間の問題を分解し、ネルダー・ミード・アルゴリズムを使用して、検索を実行する。オンラインで利用可能なＳ．Ｇ．ジョンソンによる「The Nlopt Nonlinear-Optimization Package」（http://ab-initio.mit.edu/nlopt）参照。
【００６３】
Ｃｏｂｙｌａ（Ｌ）
線形近似による位相制約下最適化（Constrained Optimization By Linear Approximation）方法は、頂点シンプレックスを用いたコスト関数および制約の線形近似を構成することによって、機能し、この近似値を最小化する。シンプレックスの半径は次第に減少し、標準形状を維持する。パウエルによる「Advances in Optimization and Numerical Analysis」（A Direct Search Optimization Method That Models the Objective and Constraint Functions by Linear Interpolation、ドルドレヒト、クルーワー アカデミック、１９９４年、５１―６７頁）、および、パウエルによる「Direct Search Algorithms For Optimization Calculations」（Dept. Appl. Math. Theoretical Phys., Cambridge Univ., U.K., Tech. Rep., 1998）参照。
【００６４】
Ｂｏｂｙｑａ（Ｌ）
二字近似による境界最適化（Bound Optimization By Quadratic Approximation）方法の各反復は、概して展補間点を考慮することによって、構成されたコスト関数の二次近似を使用する。パウエルによる「The BOBYQA Algorithm For Bound Constrained Optimization Without Derivatives」、（Centre Math. Sci., Cambridge Univ., Cambridge, U.K., Tech. Rep., 2009）参照。
【００６５】
パウエル・ブレント（Ｌ）
各反復で、一連の正確な１Ｄ線最適化は、ブレントの方法を使用して実行される。その方法は、共役（conjugate）検索方向のパウエルの方法を使用して更新される。ブレントによる「Minimization Without Derivatives, Englewood Cliffs, NJ: Prentice-Hall, 1973」、および、パウエルによる「An Efficient Method For Finding The Minimum Of A Function Of Several Variables Without Calculating Derivatives」、Comput. J.、第７巻、第２号、１５５〜１６２頁、１９６４年参照。
【００６６】
Ｐｒａｘｉｓ（Ｌ）
ブレントのＰＲｉｎｃｉｐａｌ―ＡＸＩＳ方法は、共役検索方向のパウエルの方法を更新したものであり、「Algorithms for Minimization Without Derivatives, Englewood Cliffs, NJ: Prentice-Hall, 1973」に記載されている。
【００６７】
微分展開（Ｇ）
これは、集団ベースの確率的な大域的最適化方法であり、その主要な特徴は、集団進化の基礎としての個人の対の違いのベクトルの使用である。ストーン等による「Differential Evolution-A Simple And Efficient Heuristic For Global Optimization Over Continuous Spaces」、J. Global Optimizat.、第１１巻、３４１〜３５９頁、１９９７年参照。
【００６８】
Ｄｉｒｅｃｔ（Ｇ）
この大域的最適化アルゴリズムは、本明細書の場合のように、有限境界拘束に関する課題のために設計される。パラメータ空間は、ますます小さいハイパーレクタングルに分割することによって、系統的かつ決定論的に検索される。ジョーンズ等による「Lipschitzian Optimization Without The Lipschitz Constant」、Optimizat. Theory Appl.、第７９巻、１５７〜１８１頁、１９９３年参照。
【００６９】
各アルゴリズムの自由パラメータは、アルゴリズムの原作者またはメーカーの勧告に従って調整された。以下のアルゴリズムの実施は、ＳＢＰＬＸ、Ｃｏｂｙｌａ、ＢｏｂｙｑａおよびＤｉｒｅｃｔについて、http://ab-initio.mit.edu/nloptのＮＬｏｐｔライブラリから得られる。課題の性質を変えないために、多数の境界は、オプティマイザに供給され、トランスレーションパラメータ用の±２００ｍｍ、回転パラメータ用の±４５度、スケール・パラメータ用の［0.5,1.5］、スケール・回転パラメータ用の±３６０度である[式（２）においてA1012]。
【００７０】
微分展開を除いて、基準を停止する同一セットが、全オプティマイザ用に用いられた。具体的には、オプティマイザは、１つの最適化の後、（１）コスト関数の値が、1e-14未満で減少するか、あるいは、（２）最適化されたパラメータの値が、並進パラメータに関して1e-4ｍｍ未満、回転パラメータに関して1e-5ラジアン、スケーリング・パラメータに関して1e-7だけ変更するまで、進むことが許可された。パラメータ許容度は、直観的な解釈を有する。例えば、1e-4ｍｍの並進の違いは、画像平面上の１／２０００ピクセルの最大変位に対応する。微分展開オプティマイザに関しては、確率的なアルゴリズムの収束が規則的ではないので、本発明の実施形態は、コスト関数を、約１０秒の中央値の動作に対応する１０００００回評価するまで、オプティマイザが動作するのを許容した。
【００７１】
式（１）にて表される本発明の実施形態に従うグローバル・アライメント・エネルギーを最小化する際および良質の２Ｄ／３Ｄアライメントを生成する際の９つのオプティマイザの性能は、以下の実験的な設定を使用して評価された。５人の患者のデータセットの各々に対して、３対の蛍光透視画像が考慮された。それは、ＣＴＡ獲得に時間的にアライメントされた画像および前後の隣接するフレームである。時間的アライメントは２つのモダリティ間の観察可能な違いを最小化するのを支援し、それによって、レジストレーションに失敗する危険を減らす。レジストレーション行列が毎秒ほぼ１回更新され、外科的ガイダンスを改善するため、これは臨床的に合理的である。加えて、オプティマイザのロバスト性を試験するために、１３の異なる初期点が、初期化のために使われた。T0を、装置のキャリブレーションによって与えられる変換であると仮定する。１３の初期点は、T0と、３本の主軸に沿って±５ｍｍだけ位置がずれたT0と、３本の主軸の周りに±７．５度だけ回転したT0と、である。合計５×３×１３＝１９５回の実験は、テストされる各オプティマイザに対して実行された。上述したように、式（１）は、（１）並進のみの変換、（２）剛体変換および（３）アフィン変換を連続して用いることにより、最小化可能である。
【００７２】
全９つのオプティマイザに対する３つの変換モデルの各々によって、最小化した後残った平均残余のエネルギーの平方二乗平均の分析は、アフィン変換モデルを有するエネルギー関数を最小化するための最善のアルゴリズムが以下の、微分展開、ネルダー・ミード、パウエル・ブレント、ベスト・ネイバー、ＳｂｐｌｘおよびＤｉｒｅｃｔであるということを明らかにする。平均２Ｄ射影エラーが考慮されると、順番はわずかに変化して、微分展開、パウエル・ブレント、ベスト・ネイバー、Ｓｂｐｌｘ、ネルダー・ミードおよびＤｉｒｅｃｔとなる。違いは、アライメントプロセスで使用される自動セグメンテーションと、エラーを計算するために使用される手動セグメンテーションと、の間の相違におそらく起因している。それにもかかわらず、上部にランクを付けされたアルゴリズムの性能は、満足であるように見える。ＢｏｂｙｑａおよびＣｏｂｙｌａアルゴリズムの比較的低い性能は、それぞれ、エネルギー関数の形状がそれらの二次モデルおよび線形モデルによって、良好に表されないという徴候であるかもしれない。Ｐｒａｘｉｓ方法は、この課題に最も不適切なローカル・オプティマイザであるように見える。
【００７３】
微分展開アルゴリズムが、残余エネルギーおよび平均２Ｄエラーの観点で最高の性能を有すると判明した。しかしながら、平均２Ｄエラーの増加は、ローカル・アルゴリズムより長い約１００回の計算時間に重要ではなかった。しかし、そのオプティマイザの計算時間の中央値は、ちょうど１１秒未満であり、タイム・クリティカルでないアプリケーションには適当となりうる。一般に、グローバル・オプティマイザは、最も高性能のローカル・オプティマイザより低いエネルギーまたはエラー形状に、極めてまれに至るということに留意されたい。これは、後者がグローバルな最適条件の近くにある解決法を実際に見つけるという徴候である。グローバル・オプティマイザが最低の解決法を回避することに一般に成功する点に注意する価値がある。残余エネルギーと認められた視覚の対応との間に良好な相関が存在する。それにもかかわらず、アフィン変換モデルで最高のアライメントに関してさえ、比較的大きい相違が、投影中心線と２Ｄ蛍光透視との間に存在する。これは、モダリティ間の時間的アライメントの不正確性による、データセットの非アフィン変形の存在と、さまざまな取得ステップの間、患者の位置の変化により導入される形状変化と、におそらく起因している。
【００７４】
グローバル・アライメント方法と非剛体レジストレーションとの比較
ベスト・ネイバー・オプティマイザを用いて得られた詳細なグローバル・アライメントの結果が表示され、非剛体レジストレーション後に得られた結果と比較される。全ての場合において、中心線モデルは、ＣＴ再構築用に使用される同一の心臓位相において、ゲートされる蛍光透視フレームでレジストレーションされた。並進、剛体、アフィンおよび非剛体レジストレーション後の結果は、図３および図１５の３人の患者からのデータセットに関して示され、図１５では、左から右に、入力画像、剛体アライメント、アフィンアライメントおよび非剛体レジストレーションを示す。両方の２方向画像は、すべての場合に示される。矢印は関心領域を示す。２Ｄマニュアル・セグメンテーションに関する２つの透視面上の中心線の平均残余の２Ｄ射影エラーも、５人の患者に対する剛体、アフィンおよび非剛体レジストレーション後に、（ミリメータで）算出され、その結果は図１６の表Ｉに示されている。図１５の上の行および下の行は、それぞれ、表Ｉの患者１および患者２に対応する。図１５の２つの左下の画像フレームのＣＴＯラベルは、慢性完全閉塞を意味する。非剛体レジストレーションに必要な計算時間は、７３０ｍｓ〜３３００ｍｓの範囲にあり、平均１５９１ｍｓであった。
【００７５】
モデルの複雑さが増加するにつれて、平均２Ｄ射影エラーおよびそれらの標準偏差は減少するか、または、ほぼ一定のままである。しかしながら、剛体、アフィンおよび非剛体変換モデルの相対的な貢献は、データセットごとに変化する。これは、個々のケースにより示される変形の性質によって説明可能であり、変形の性質は、患者の位置、ＣＴＡと２Ｄ蛍光透視獲得との間の間隔および獲得プロトコルに関連している。また、患者の呼吸によって、非剛体の心臓の変形が生じうる。心臓の鼓動は非剛体変形の原因でもあるため、ＥＣＧゲーティングによる時間的アライメントの精度も重要である。最大エラーもモデル複雑さに伴って減少する傾向があるが、時々安定なままである。このような場合、自動のセグメンテーション方法により計算されるため、レジストレーション・アルゴリズムの性能は不完全な２Ｄセグメンテーションにより制限されるかもしれない。実際、ほとんど情報が利用できない領域では、本発明の実施形態による非剛体レジストレーション方法は、中心線をできるだけ変えてはならない。この特性は、慢性完全閉塞（ＣＴＯ）の場合のように、構造の一部がＸ線蛍光透視画像上で十分には見えない状況において望ましい。例えば、図１５の下行は、中心線がＣＴＯの上下に良好にアライメントされているが、コントラストが低い部分にほとんど変形がないことを示す。
【００７６】
マルチ・フレーム・シナリオのグローバル・アライメント
２方向の一対に対する中心線のアライメントの他の潜在的に役立つシナリオは、２方向血管造影法からの時間的フレームシーケンスに対する中心線のアライメントである。マルチ・フレーム・シナリオにおける本発明の一実施形態によるグローバル・アライメント方法の性能は、患者１のＬＣＡデータセットから一連の１１の２方向フレームを使用して示される。グローバル・アライメントは、本発明の実施形態に従う剛体およびアフィン変換モデルを使用して計算された。式（５）のフレーム間組織化パラメータγは、０．０〜１．０の間で変化した。式（５）のコスト関数を最小化するのに最も効果的であると判明したので、ベスト・ネイバー・オプティマイザがこの場合使われた。得られたサンプル結果は図１７に示され、図１７はフレームシーケンスの上のレジストレーションを示す。図１７において、一番上の行は、γ＝0.1およびEMulti*＝120.284を有する剛体変換モデルを表し、中間の行は、γ＝0およびEMulti*＝119.356を有するアフィン変換モデルを表し、一番下の行は、γ＝0.1およびEMulti*＝116.905を有する正規化アフィン変換モデルを表し、ここで、EMulti*は、式（５）の右辺の第１項に対応する。矢印は関心領域を示す。計算時間は、一番上の行、中間の行および一番下の行で示される実験に対して、それぞれ、３９秒、８９秒および９１秒であった。予想されるように、アフィン変換モデルを使用すると剛体変換モデルを使用するより良好な定性的結果に至ることが判明した。アフィン変換モデルを使用するとき、組織化（γ＝０．１）を使用することも結果を改善した。大きい組織化値（γ＞１．０）を用いると、モデルの剛性は増加し、オプティマイザは初期位置から遠い極小値を見つけることができない。予備調査は、例えば、パラメータが協調して変化可能なことによって、領域に特有のオプティマイザがより良好な結果を達成しうることを示唆する。それにもかかわらず、本発明の実施形態に従う設定は、大部分の心周期の間、ＬＣＡのアライメント変換をトラッキングするのに適切に聞こえた。ＲＣＡデータセットに関する実験は、説得力のある結果につながらなかった。事実、本発明の一実施形態によるグローバル・アライメント手順が機能しても、ＣＴＡおよびＸ線取得が近い時間的アライメントに近い場合を除き、心周期の間、ＲＣＡにより支持される剛体変形の量によって、アフィン・モデルの効果がなくなる。
【００７７】
右冠動脈の半自動のトラッキング
この実験では、本発明の実施形態による非剛体レジストレーション方法に基づく半自動手順を用いて、ＲＣＡのトラッキングを行った。ゲートされたフレームから開始して、本発明の一実施形態によるグローバル・アライメントおよび非剛体レジストレーションが実行される。次に、本発明の一実施形態によるこの変形した中心線モデルは、次の一対の蛍光透視フレームのための第一のモデルおよび位置としてオペレータにより用いられる。図１８に示されるように、このプロセスは１つの心周期上の全フレームのために繰り返され、図１８は、左から右に、ゲートされたフレーム、＋４、＋６、＋８および＋１０フレームを表す。第一の曲線は１８１であり、レジストレーションされた曲線は１８２である。了解されるように、この手順は、１つの心周期の間、ＲＣＡのトラッキングを成功させる。この実験がシミュレーション研究の前に行われたので、わずかに異なるパラメータセット、μ＝０．０５、υ＝５．０およびλ＝０．１が用いられた。必要な時間は約１５分であり、これはオフラインの設定において、合理的である。
【００７８】
考察
臨床データの実験に基づいて、ネルダー・ミード、パウエル・ブレント、ベスト・ネイバーおよびＳｂｐｌｘのローカル・オプティマイザが短時間（中央値＜１０５ｍｓ、最大値＜９５０ｍｓ）で良好な結果を与える。シミュレーションも考慮すると、ネルダー・ミード・アルゴリズムは、全体的に最善の性能を与えた。２つのグローバル・オプティマイザは、その結果が大きな改良につながることは極めてまれであり、一桁高い計算時間を必要とした。計算時間を無視すると、微分展開アルゴリズムは、一般に、最善の解決案を提供した。ローカル・オプティマイザを使用するとき、アライメント時間は一貫して１秒未満であるため、この方法は、インターベンション中の使用に適するものになる。剛体変換の代わりにアフィン変換を使用することの利点は、データセットの性質に依存し、場合によっては重要になりうる。３Ｄおよび２ＤモダリティはＥＣＧゲーティングを使用して時間的にアライメントされるとき、全体として、グローバル・アライメント手順がＬＣＡおよびＲＣＡデータセットの両方の使用に適切であること判明した。ＣＴＯ患者からのデータセットに関する実験は、類似の臨床ケースに適用されるとき、本発明の実施形態に従う方法の利点を示す。本発明の実施形態による方法は、ＣＴＡ獲得からセグメント化された３Ｄ形状を有する２Ｄ画像を増加することによって、インターベンション・ガイダンスを支援することもでき、このことにより２Ｄ画像の解釈に固有の曖昧性を著しく低減する。
【００７９】
非剛体レジストレーションに関しては、本発明の一実施例による変形モデルに含まれる正規化は、中心線を最小量だけ変形させることによって、残りの構造との一貫性を確実に保つことによって、欠測データ（例えば、図１５の患者２）を有する領域の管理を可能にする。過剰な変形または矯正は発生しない。全体の計算時間は、通常、３秒未満であり、これはインタラクティブなアプリケーションに使用できる。本発明の一実施形態による非剛体レジストレーション方法は、適切な初期点から開始するとき、隣接する２Ｄ構造に対する中心線をスナップし、特に困難なローコントラストＣＴＯケースに関して、２Ｄ画像の解釈をより容易にする。これも、操作上の画像と一緒に更新された３Ｄモデルを示すことを可能にし、３Ｄ空間の改良された認識を提供できる。
【００８０】
最終的に、複数のＸ線２方向血管造影法フレームを使用する実験によって、本発明の一実施形態によるマルチフレーム・グローバル・アライメント方法が、ＬＣＡデータセットに良好に作用するとともに、本発明の一実施形態によるフレーム間の正規化が、改良された結果につながるということが判明した。本発明の一実施形態による非剛体半自動のトラッキング方法は、より大きな非剛体変形を有するケース、例えば非時間的にアライメントされたＲＣＡデータセットを有する場合を扱い、そのようなデータセットに適用された。本発明の一実施形態による半自動の方法は、１つの心拍をカバーしているシーケンスに対して約１５分の相互作用を必要とし、それは実用的なオフラインのアプリケーションに合理的に思える。
【００８１】
システムの実現
本発明の実施形態は、各種形態のハードウェア、ソフトウェア、ファームウェア、特殊目的のプロセスあるいはその組み合わせによって実現可能であることを理解されたい。一実施形態では、本発明は、コンピュータ可読のプログラム記憶装置に具体的に組み込まれたアプリケーションプログラムとしてソフトウェアで実現可能である。アプリケーションプログラムは、任意の適切な構造を有する機械にアップロード可能であるとともに、当該機械によって実行可能である。
【００８２】
図１９は、ライブ蛍光透視画像を用いた冠動脈モデルの非剛体２Ｄ／３Ｄレジストレーションを実現するコンピュータシステムのブロック図である。図１９を参照すると、本発明を実施するためのコンピュータシステム191は、とりわけ、中央処理装置（ＣＰＵ）192と、メモリ193と、入出力（Ｉ／Ｏ）インタフェース194と、を具えることができる。コンピュータシステム191は、一般的に、Ｉ／Ｏインタフェース194を介してディスプレイ195と、マウスやキーボードのような複数の入力装置196と、に接続されている。支持回路は、キャッシュ、電源、クロック回路および通信バスのような回路を含むことができる。メモリ193は、ランダムアクセスメモリ（ＲＡＭ）、読み出し専用メモリ（ＲＯＭ）、ディスクドライブ、テープドライブ等、あるいはその組み合わせを含むことができる。本発明は、ルーティン197として実現され、このルーティン197はメモリ193に記憶され、ＣＰＵ192によって実行され、信号源198からの信号を処理する。それゆえ、コンピュータシステム191は、汎用コンピュータシステムであり、本発明のルーティン197を実行するとき特殊目的のコンピュータシステムになる。
【００８３】
コンピュータシステム191は、また、オペレーティングシステムと、マイクロ命令コードと、を含む。本明細書に記載した様々な処理および機能は、オペレーティングシステムを介して実行されるマイクロ命令コードの一部あるいはアプリケーションプログラムの一部（あるいはその組み合わせ）とすることができる。さらに、様々なその他の周辺装置、例えば、追加のデータ記憶装置および印刷装置を、コンピュータプラットフォームに接続することができる。
【００８４】
添付図面に示されたシステムの構成要素および方法のステップのいくつかは、ソフトウェアによって実施可能であるため、システムの構成要素（あるいはプロセスステップ）間の実際の接続は、本発明がプログラムされる方法に依存して異なるということを理解されたい。本明細書に記載した本発明の教示を考慮すると、当業者は本発明と同一および類似の実施形態あるいは構成を考えることができる。
【００８５】
本発明は、例示的実施形態を参照して詳細に説明してきたが、当業者であれば、添付の特許請求の範囲に規定した本発明の精神および範囲から逸脱することなく各種変更および置換が可能であることを理解するであろう。
【技術分野】
【０００１】
米国出願の相互参照
この出願は、２０１１年９月２８日にハリ・サンダーにより出願された米国仮出願第６１/５４０，１３４号の「心臓インターベンションのガイダンス用のライブ蛍光透視法を用いた冠動脈モデルの非剛体２Ｄ／３Ｄレジストレーション」の優先権を主張するものである。
【０００２】
本発明は、デジタル医用画像の解剖モデルの非剛体２Ｄ／３Ｄレジストレーションの方法に関するものである。
【背景技術】
【０００３】
Ｘ線透視法は、多数の低侵襲的処置、例えば、慢性完全閉塞（ＣＴＯ）の経皮冠動脈インターベンション（ＰＣＩ）のガイダンス用に選択されるモダリティである。これらの処置の間、ガイドワイヤを用いたＣＴＯの横断（crossing）は特に有害である。なぜなら、閉塞が造影剤の伝播を阻止し、蛍光透視法下において、血管の閉塞部を不可視にするからである。さらに、蛍光透視法の射影性質は、３Ｄ構造の解釈を不明確にし、そのことがさらにインターベンションを複雑にする。同様の問題は、低侵襲的処置にも無数に存在している。
【０００４】
現在のプレインターベンショナル計画は、通常、ＣＴ血管造影（ＣＴＡ）あるいは他の３Ｄ画像モダリティの取得を含み、視覚的な不確実性を低減している。血液プール造影剤（contrast）注入を用いたＣＴＡの性能特性は、直接の造影剤注入を用いたインターベンショナルＸ線透視法によって認識可能なものと対照的に、ＣＴＯを引き起こす石灰化を明らかに区別可能にする。これは、計画段階にてＣＴＡを非常に有益なツールとするが、これらのデータとインターベンショナル画像との間の対応を確立することは挑戦すべき問題であると判明している。
【０００５】
この問題に対処するために、処理前の３次元モデルは取得されたボリュームから抽出され、２Ｄ蛍光透視図とアライメントされ、ライブ画像上に重ねられ、それによって、インターベンショナル画像を増加させることができる。しかしながら、このアライメントを達成することは本来困難である。なぜなら、主には、インターモダル対応を見つけることは些細ではないためであり、根本的な最適化の非線形態様のためでもある。さらに、満足な２Ｄ／３Ｄレジストレーションを提供するために、単純な剛体変換では十分でないおそれがある。３Ｄ計画画像が呼吸停止（breathold）の状態で得られるので、自由な呼吸状態で得られる術中の画像と比較して、著しい形状変化が存在する。これは、術前モデルを術中のライブ画像にアライメントする間、非剛体レジストレーション方法を使用することを役立つようにする。しかしながら、これは、困難な設定不備の逆の作業である。それにもかかわらず、２Ｄ／３Ｄレジストレーション方法には、既存の臨床の流れに最小の変更を加えるだけで、インターベンショナルＸ線血管造影法に関する不確定性を大きく低減する可能性がある。より広い展望から、２Ｄ／３Ｄレジストレーション方法は、例えば神経学、整形外科および心臓学の分野における多数の用途を有する。
【０００６】
２Ｄ／３Ｄレジストレーションの最適化面は、画像の離散化、レジストレーションされる構造の複雑さおよび剛体あるいはアフィン変換モデルの使用といったいくつかの理由のために非常に非線形になりうる。その結果、コスト関数の最小化は、困難である。レジストレーションプロセスによって、潜在的に達成可能な最大精度は、含まれる変換モデルの複雑さに直接関連している。多数の２Ｄ／３Ｄレジストレーション方法は、剛体変換モデルあるいはわずかにフレキシブルなアフィン・モデルのみを考慮している。一般的に、これは、骨のような剛体構造や、モダリティの第一のアライメントを与えるには適切であるが、フレキシブルな構造の形状変化を説明するには不十分であると判明している。したがって、非剛体変形モデルは、２Ｄ／３Ｄレジストレーションに対して提案されてきた。
【０００７】
臨床状況において、１つまたは２つの透視面は、インターベンション・ガイダンスのために用いられる。モノプレン（Monoplane）のシナリオが研究されてきたが、報告された誤差は、面外方向において大きい。この正確性の欠如は、レジストレーションプロセスにおける信頼性を著しく制限する。２方向の取得は、これらのインターベンションに関連した不明確性を大きく低減し、今日、日常的に複雑な心臓インターベンションに実行されている。
【発明の概要】
【０００８】
本願明細書に記載されている本発明の実施形態には、概して、２Ｄ／３Ｄ非剛体レジストレーションためのシステムと方法であって、術前の冠動脈の３Ｄ中心線モデルを術中の２以上の蛍光透視画像に非剛体でアライメントし、インターベンショナル画像を視覚的に増加させるシステムと方法が含まれる。本発明の実施形態によるレジストレーションされたモデルは、インターベンショナル血管造影図の上に重ねられ（オーバーレイ可能）、画像導入の慢性完全閉閉塞処理の間、外科的援助を提供し、２Ｄインターベンショナル画像に固有な不確定性を低減する。本発明の一実施形態による方法は、２つの部分、すなわち、（１）第一の固体またはアフィンアライメントを提供するために算出されるグローバル変換モデルと、（２）最終的なレジストレーションを計算するために使用される完全非剛体モデルと、を含む。グローバル・アライメント・パラメータを推定するために、距離マップを用いて、３次元モデルの射影と２Ｄ画像から自動抽出される特徴との間の相違を測定する。第二の部分では、完全非剛体レジストレーション方法を用いて、任意の局所的な形状相違を補償する。本発明の一実施形態による非剛体レジストレーション方法は、変化のフレームワークに基づき、同時照合および再構築プロセスを使用して、非剛体レジストレーションを計算する。最新のハードウェア上の３秒未満の一般的なランタイムでは、本発明の一実施形態によるアルゴリズムは、ＣＴＯ手順の間、手術中に使われるのに十分急速に解析可能である。さらに、本発明の一実施形態による２Ｄ／３Ｄレジストレーション方法を用いて、フレームシーケンス全体にわたって、対象構造の動きをキャプチャすることもできる。
【課題を解決するための手段】
【０００９】
本発明の一実施形態によれば、心臓インターベンション中の２Ｄ蛍光透視画像を用いたデジタル３Ｄ冠動脈モデルの非剛体レジストレーションのための方法が提供され、前記方法は、Qsの３Ｄ制御点からなる一組のSセグメントを具える３Ｄ中心線であって、冠動脈ツリーのデジタル化された３Ｄ中心線の表示を提供するステップと、前記３Ｄ中心線を、少なくとも２つの２Ｄ蛍光透視画像に、グローバル・アライメントするステップと、エネルギー汎関数を最小化することによって、前記３Ｄ中心線を、少なくとも２つの２Ｄ蛍光透視画像に、非剛体レジストレーションするステップと、を含み、前記エネルギー汎関数は、再構築された中心線点とレジストレーションされた中心線点との間の二乗差の総和と、二乗３Ｄレジストレーション・ベクトルの総和と、二乗微分係数３Ｄレジストレーション・ベクトルの総和と、心筋分岐エネルギーと、を含み、前記３Ｄ中心線の非剛体レジストレーションは、前記３Ｄ中心線の座標系にて、対応する中心線点xs,qに適用される一組の３Ｄ並進ベクトルrs,qとして表される。
【００１０】
本発明のさらに別の一実施形態では、前記冠動脈ツリーは前記２Ｄ蛍光透視画像からセグメントされ、前記セグメントされた冠動脈ツリーは前記２Ｄ蛍光透視画像の各々の２値化として表される。
【００１１】
本発明のさらに別の一実施形態では、前記３Ｄ中心線を、前記少なくとも２つの２Ｄ蛍光透視画像に、グローバル・アライメントするステップは、前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線の並進（translation）のみのレジストレーションを実行するステップと、前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線の剛体レジストレーションを実行するステップと、前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線のアフィン・レジストレーションを実行するステップと、
を含み、各ステップでは、
【数１】
の形のエネルギー汎関数は最小化され、ｊは、前記並進のみのレジストレーションと、前記剛体レジストレーションと、前記アフィン・レジストレーションと、を分類し、Cは、前記３Ｄ中心線の座標系と３Ｄ画像形成装置の座標系との間の座標変換であり、Tjは、前記並進のみのレジストレーションと、前記剛体レジストレーションと、前記アフィン・レジストレーションと、に対して、前記３Ｄ画像形成装置の座標系と、前記蛍光透視画像を取得する２方向Ｃアームの座標系と、の間の座標変換であり、xは、均一な座標の３Ｄ点であり、nは、前記蛍光透視画像上の合計を表し、sおよびqは、前記３Ｄ中心線の前記セグメント上の制御点を分類し、Pnは、前記２方向Ｃアームの前記座標系と、第ｎ番目の蛍光透視画像の前記座標系と、の間の座標変換であり、Ψnは、前記３Ｄ位置ｘを前記第ｎ番目の蛍光透視画像に射影する射影オペレータであり、
【数２】
は、前記第ｎ番目の蛍光透視画像の２Ｄ点と、前記第ｎ番目の二値化された画像Bnの前記セグメントされた冠動脈ツリーに最も近い点と、の間の距離である。
【００１２】
本発明のさらに別の一実施形態では、変換行列Tjは、再センタリングされ、前記変換行列Tjの回転パラメータの変化は、前記３Ｄ画像形成装置の前記座標系の原点とは対照的に、前記３Ｄ中心線の前記座標系の原点の周りで発生するように見える運動を誘導する。
【００１３】
本発明のさらに別の一実施形態では、時間的に連続したF対の２Ｄ蛍光透視画像を提供するステップと、前記エネルギー汎関数
【数３】
を最小化することによって、前記３Ｄ中心線を、各対の２Ｄ蛍光透視画像にグローバル・アライメントするステップと、をさらに含み、fは時点の合計であり、Tfは各時点の変換であり、γは所定のパラメータであり、D(Tf・C・xs,q,Tf+1・C・xs,q)はフレーム間距離である。
【００１４】
本発明のさらに別の一実施形態では、前記３Ｄ中心線を少なくとも２つの２Ｄ蛍光透視画像にグローバル・アライメントするステップの後、前記３Ｄ中心線の点を、前記２Ｄ蛍光透視画像の対応する血管中心線点に照合し、前記照合した２Ｄ点から３Ｄ中心線点
【数４】
を再構築するステップをさらに含む。
【００１５】
本発明のさらに別の一実施形態では、前記点照合は、各３Ｄ中心線点xs,pの射影x'を、前記グローバル・アライメントからの各２Ｄ蛍光透視画像上に見つけるステップと、x'から開始して、各２Ｄ蛍光透視画像に対応するセグメンテーション画像Bnを、前記セグメンテーション画像B上の前記３Ｄ中心線の射影に垂直な線に沿って検索し、中心線点xs,pに照合する２つまでの点を見つけるステップと、xs,pに最も近い点を照合点xnとして選択するステップと、前記蛍光透視画像の照合点を用いて、前記心臓インターベンションに用いられている画像処理システムの射影形状を使用している３Ｄ中心線点
【数５】
を再構築するステップと、を含む。
【００１６】
本発明のさらに別の一実施形態では、再構築された中心線点
【数６】
とレジストレーションされた中心線点
【数７】
との二乗差の前記総和は、
【数８】
であり、
【数９】
は指標関数
【数１０】
であり、照合点がいかなる画像平面でも見つからない場合、あるいは、２Ｄ蛍光透視画像上の
【数１１】
の射影x'と前記第ｎ番目の２Ｄ蛍光透視画像の照合点xn"との間の２Ｄ距離が所定の閾値Max2Dより大きい場合、あるいは、中心線点xs,qと対応する再構築された点
【数１２】
との間の３Ｄ距離が所定の閾値Max3Dより大きい場合、H(xs,p)=0であり、その他の場合、H(xs,p)=1である。
【００１７】
本発明のさらに別の一実施形態では、二乗３Ｄレジストレーション・ベクトルの前記総和と、二乗微分係数３Ｄレジストレーション・ベクトルの前記総和と、前記エネルギー汎関数の前記心筋分岐エネルギーは、
【数１３】
であり、rは、レジストレーション・ベクトルであり、xは、前記３Ｄ中心線上の制御点であり、sおよびqの前記総和は、前記３Ｄ中心線の全セグメントの全制御点の合計であり、Kは、前記３Ｄ中心線の接合部から距離Ds離れた各セグメントの点を、前記接合部で接合される他のセグメントの対応する制御点に連結することにより形成される制御点の全対の一組であり、前記接合部では、連結がレジストレーションの間拡大または圧縮される場合、各連結はその変位に比例した力を発生させる。
【００１８】
本発明のさらに別の一実施形態では、前記エネルギー汎関数は、
【数１４】
を反復することによって、最小化され、zは、点rs,p=[0,0,0]T|∀{s,q}から反復z=0で開始する反復カウンタであり、εは、小さい定数であり、Ks,q⊂Kは、制御点(s,q)に作用する一組の制約であり、μ、υおよびλは、所定のパラメータである。
【００１９】
本発明のさらに別の一実施形態では、コンピュータ可読の非一時的プログラム記憶装置が提供され、前記コンピュータにより実行されるとともに、心臓インターベンション中の２Ｄ蛍光透視画像を用いたデジタル３Ｄ冠動脈モデルの非剛体レジストレーションのための方法を実行する命令プログラムを有する。
【図面の簡単な説明】
【００２０】
【図１】本発明の一実施形態に係る、対象とするシステムの形状の概略図である。
【図２】（ａ）〜（ｈ）は、発明の一実施形態に係る、レジストレーションの複数のステージによる入力データの変換を示す図である。
【図３】発明の一実施形態に係る、非剛体レジストレーションの進行を示す図である。
【図４】発明の一実施形態に係る、中心線制御点照合アルゴリズムのフローチャートである。
【図５】（ａ）（ｂ）は、本発明の一実施形態に係る、点照合および分岐点における分岐制約の形成を示す。
【図６】本発明の一実施形態に係る、反復回数の関数として、中心線の全移動を示すグラフである。
【図７】本発明の一実施形態に係る非剛体レジストレーション方法において、パラメータλの挙動を示す図である。
【図８】本発明の一実施形態に係る非剛体レジストレーション法の後、平均２Ｄおよび３Ｄ誤差を示す図である。
【図９】本発明の一実施形態に係る、非剛体２Ｄ／３Ｄレジストレーション法のフローチャートである。
【図１０】（ａ）〜（ｄ）は、本発明の一実施形態に係るＤＤＲの形成を示す図である。
【図１１】本発明の一実施形態に係る初期位置の摂動に関して、オプティマイザの性能を例示する。
【図１２】本発明の実施形態に係るノイズ・レベルに対する入力ＤＲＲノイズの標準偏差に関して、オプティマイザの性能を例示する。
【図１３】本発明の実施形態に係るシミュレーションデータを有する若干のサンプル非剛体レジストレーションを示す。
【図１４】本発明の実施形態に係るシミュレーションの非剛体変形レベルに関する３Ｄ残差エラーのグラフである。
【図１５】本発明の実施形態に係る、剛体、アフィンおよび非剛体レジストレーション後の入力画像および結果を表す。
【図１６】本発明の実施形態に係る、５人の患者のための剛体、アフィンおよび非剛体レジストレーション後の２Ｄマニュアル・セグメンテーションに関して２つの透視面上の中心線の平均残余の２Ｄ射影エラーを示している表である。
【図１７】本発明の実施形態に係る、一連のフレーム上のレジストレーションを示す。
【図１８】本発明の実施形態に係る、１つの心周期の複数のフレームによるＲＣＡのトラッキングを表す。
【図１９】本発明の一実施形態に係る、ライブ蛍光透視画像を用いた冠動脈モデルの非剛体２Ｄ／３Ｄレジストレーションのための方法を実施するためのコンピュータシステムの一例のブロック図である。
【発明を実施するための形態】
【００２１】
本発明の例示的実施形態は、概して、ライブ蛍光透視画像を用いた冠動脈モデルの非剛体２Ｄ／３Ｄレジストレーションのためのシステムおよび方法を含む。したがって、本発明は各種変更および代替実施形態の影響を受けるが、特定の実施形態を、単に一例として図面および明細書に詳細に記載する。しかしながら、本発明は、記載されたこの特定の形態に限定することを意図するものではなく、本発明は、本発明の精神および範囲内にある全ての変更、等価および代替を含むものである。
【００２２】
本明細書に記載の用語「画像」は、個別の画像要素（例えば、２次元画像のピクセルおよび３次元画像のボクセル）から構成された多次元データを意味する。画像は、例えば、コンピュータ断層撮影法、磁気共鳴映像法、超音波、あるいは、従来当業者に知られている他の任意の医用画像システムによって収集された被写体（subject）の医用画像とすることができる。画像はまた、医学とは関係ない分野、例えば、リモートセンシングシステム、電子顕微鏡法等から提供される場合もある。画像は、R3からRまたはR7の関数とみなすことができるが、本発明の方法はこれらの画像に限定されず、任意の次元、たとえば、２次元画像あるいは３次元体積の画像に適用可能である。２次元あるいは３次元画像のために、画像の定義域は一般的に２次元あるいは３次元の長方形配列であり、各ピクセルあるいはボクセルは、２本あるいは３本の相互に直交する軸に関して方向付け可能である。本明細書に記載の用語「デジタル」および「デジタル化」は、デジタル取得システムを介して、あるいはアナログ画像からの変換によってデジタル化されたフォーマットにおける画像あるいは体積を適切に意味する。
【００２３】
本発明の一実施形態によると、術前のＣＴＡからセグメントされた、冠動脈ツリーの３Ｄ中心線表現は、２以上の同時の蛍光透視画像に非剛体レジストレーション可能である。初期設定の３Ｄ位置から開始し、装置のキャリブレーションから計算し、中心線の３Ｄアライメントは、並進のみの運動と剛体変換とアフィン変換とを用いて徐々に調節される。最後に、非剛体変換が計算され、最終的なレジストレーションが提供される。
【００２４】
本発明の一実施形態によると、冠動脈ツリーは、制御点xs,p∈R、s∈{1,...,S}、q∈{1,...,QS}から構成された一組のSセグメントQsによって表され、方向付けのない（undirected）非巡回グラフを形成する。説明を明確にし、一般性を失わないために、以下の説明は、１対の蛍光透視画像In、n∈{1,2}が与えられる場合に制限されている。しかしながら、この制限は例示であり限定的ではなく、本発明の実施形態に係る方法は、複数の蛍光透視画像に拡張可能であることを理解されたい。さらに、説明を明確にし、一般性を失わないために、両者が同一のサイズW×H∈N2であると仮定する。血管は、既知の方法を用いて入力蛍光透視画像から自動的にセグメント化され、入力画像In(i,j)の２値化Bn(i,j)として表される。位置(i,j)∈{0,...,W-1}×{0,...,H-1}における入力画像ピクセルが対象の構造に対応する場合、Bn(i,j)=1であり、その他の場合、Bn(i,j)=0である。
【００２５】
図１には、対象とするシステムの形状が概略的に示され、５つの座標系（ＣＳ）を用いて記載されている。図１中のラベルつきのフレームは座標系を表し、破線および太字の符号は変換行列を表す。冠動脈ツリーを表す３Ｄ中心線はＵに関して記載され、座標系（ＣＳ）Ｃは３Ｄ画像装置の基準点に中心がある。２Ｄ血管造影法に関して、Ｘは２方向Ｃアームの基準ＣＳであり、ＣＳのP1およびP2は、２Ｄ結像面に中心を置いている。ＣとＵとの間の変換およびＸと{P1,P2}との間の変換は、装置のキャリブレーションから知られており、それぞれ、剛体変換行列C、P1およびP2で符号化される。２つのＸ線平面の射影形状もまた利用可能に仮定され、射影オペレータΨ1およびΨ2で符号化される。射影オペレータΨ1およびΨ2は、同次座標x=[x,y,z,1]Tの３Ｄ点を、対応する透視面上の２Ｄ点x’=[x,y]にマップする。２つの画像診断座標系の間の変換は、Tによって表される。
【００２６】
完全性のために、例示であり限定的ではないが、自動セグメンテーション方法を以下に簡単に説明する。（１）入力蛍光透視画像は、ヘシアンベースのサトー血管フィルタ（Hessian-based Sato vesselness filter）を用いて処理され、（２）血管測定の閾値化をまばらにサンプリングすることによって、シードポイントが生成され、（３）これらのポイントから開始し、各ピクセル位置にて計算されたヘシアン行列の第１の固有ベクトルに沿って統合されることによって、ファイバが生成され、（４）ファイババンドルは反復して間引かれ（thinned）、血管の中心線を引き出す。
【００２７】
並進、剛体およびアフィン変換
一般的な臨床状況では、３Ｄ冠動脈中心線をＣＴＡ取得から２方向透視における同一構造の画像にアライメントする変換は、大まかにしかできない。なぜなら、患者の位置は高い精度で制御することができないためである。本発明の一実施形態に係る変換T、すなわち、２つの画像診断座標系間のマッピングは、装置のキャリブレーションによって推定可能であり、対象の構造の形状を考慮することによって調整可能である。初期条件から開始し、改良されたアライメント変換（alignment transformations）は、最小化プロセスを用いて計算され、グローバル変換モデルあるいは剛体変換モデルのパラメータが推定される。この点から、アフィン変換モデルもまた用いられ、３Ｄ中心線を変形し、ＣＴＡと２方向Ｘ線との間の形状の不一致を補償する。
【００２８】
本発明の一実施形態によると、３Ｄ中心線の射影と、各透視面からセグメント化された血管中心線との間の総距離は最小化される。χ(x)=T・C・xは、Ｕに関する点x=[x,y,x,1]をX座標系にマップするアフィン変換演算子とする。また、Пn(x)=Ψn・Pn・xは、Xからの点を透視面Pnにマップする射影である。本発明の実施形態は、エネルギーを、２つの蛍光透視画像で最小化され加算されるべき量として、以下のように定義する。
【数１５】
ここで、
【数１６】
は２Ｄ点とBn(i,j)=1である最接近点との間の距離である。各画像位置に対する距離
【数１７】
は、各Bn(i,j)セグメンテーションの距離変換を形成することによって、事前に計算される。２つの蛍光透視画像を使用することは例示であり限定ではなく、２つより多い蛍光透視画像が用いられる場合、加算は容易に拡張されることを理解されたい。
【００２９】
図２（ａ）〜（ｈ）は、入力データの変換を示している。図２（ａ）はＣＴ体積からの湾曲した（curved）スライスを示し、図２（ｂ）は３Ｄ中心線を示し、図２（ｃ）（ｄ）は、異なる視点からの透視血管造影図を示し、図２（ｅ）（ｆ）は蛍光透視画像のセグメンテーションを示し、図２（ｇ）（ｈ）はセグメンテーションの距離変換を示す。
【００３０】
Tの１２の値を直接最適化することは、効果的でなく非効率なため、本発明の実施形態では、代わりに、パラメトリック手法を用いる。並進のみの変換TT、剛体変換TRおよびアフィン変換TAは、それぞれ、パラメータセットT∈R3、R∈R6およびA∈R12を用いて表される。TSのパラメトリック表現から行列表現へのマッピングは、以下のように定義される。
【数１８】
【数１９】
【数２０】
【００３１】
ここで、Sijは、Sの添え字{i, ... ,j}を伴う要素を含む部分集合である。
【数２１】
【数２２】
【数２３】
ここで、cx=cos(θx)、sx=sin(θx)等である。本発明の実施形態では、変換行列TT、TRおよびTAは、再センタリングされる。回転パラメータのいかなる変化も、Ｃの原点と対照的に、Ｕ座標系の原点周辺で発生するように見える運動を誘導するというような方法で、この動作は、式（１）のエネルギー関数に影響しないが、通常オプティマイザに有益である。なぜならこれは、直観的に、ＣＴボリュームの重心で発生するように見える回転動作またはスケーリング動作に関連して、この種の再センタリングが、パラメータ空間において、距離を減らすためである。最終的なTの定義は、式（４）のとおりである。
【数２４】
ここで、ＣはＵからＣへの変換であり、j∈{T,R,A}であり、選択された変換モデルに依存する。式（１）は、３つの変換モデルの各々に対して、式（１ａ）の形をとなる。
【数２５】
任意の非派生的な多目的最適化アルゴリズムを用いて、式（１）を最小化できる。
【００３２】
初期設定の患者の位置から開始して、本発明の実施形態は、（１）並進のみの変換、（２）剛体変換および（３）アフィン変換を考慮することによって、アライメントを徐々に調整する。
【００３３】
マルチ・フレーム・アラインメント
本発明の一実施形態に係るグローバル・アライメント方法は、マルチ・フレーム・シナリオをカバーするために拡大されることができる。この場合、実施形態は、時間的連続性のため、推定されたアライメント変換が１フレームからその他にほとんど何も変えないと仮定する。本発明の一実施形態に係るグローバルなエネルギー関数は、一連の２方向フレームのために以下のように、再定義可能である。
【数２６】
あるいは
【数２７】
である。ここで、γは自由パラメータであり、χfは、χと同様に定められるが、時点ごとの異なる変換行列Tfを有する。変換行列Tfは、上述したように、並進のみ、リジッドあるいはアフィンとすることができる。上述したように、２つの蛍光透視画像上の合計nは例示であり限定的ではなく、２つ以上の蛍光透視画像が使われている場合、合計は実施形態において、容易に拡張されうる。
【００３４】
フレーム間距離D(χf,χf+1)は、リジッド変換あるいはアフィン変換とすることもできる。リジッド変換行列およびアフィン変換行列の空間が線形でないので、本発明の一実施形態に係るフレーム間距離D(χf,χf+1)の定義は精巧である。本発明の実施形態により仮定されているように、剛体変換空間は方向における小さい変化のために局所的に線形化可能である。このフレームを用いて、フレームfからフレームf+1への剛体変換は、並進ベクトル
【数２８】
および回転ベクトル
【数２９】
によって表され、
【数３０】
となるように、単位ベクトル
【数３１】
および回転角度θf,f+1のプロダクトとして定義される。並進ベクトルは、剛体変換行列
【数３２】
の並進部分に対応し、行列
【数３３】
と回転ベクトル
【数３４】
との間のマッピングは、
【数３５】
および
【数３６】
によって与えられる。
【００３５】
この表現を使用して、実施形態は、以下の通りにフレーム間距離を形式的に定めることができる。
【数３７】
【数３８】
ここで、αは、回転部分および並進部分のそれぞれの貢献のバランスをとる実数である。本発明の実施形態は、その内容全体が本願明細書に参照によって取り込まれる、ボワヴェール（Boisvert）等による「Geometric Variability Of The Scoliotic Spine Using Statistics On Articulated Shape Models」（ＩＥＥＥ Ｔｒａｎｓ. Ｍｅｄ. Ｉｍａｇ.、第２７巻、第４号、５５７〜５６８頁、２００８年４月）に従って、αを０．０５に設定する。
【００３６】
アフィン変換行列の間の線形距離を定めることは、関連が深い複雑な計算を必要とする。本発明の実施形態では、剛体変換が変換の大部分を捕えると仮定されるので、本発明の実施形態は、変換の剛体部分を正規化し、残りのアフィン・パラメータを無拘束のままにする。形式的に、アフィン変換の場合、正規化定義の変換TfはパラメータA16を使用して構成されるが、１２のパラメータのフルセットはエネルギー定義の変換χfを構成するために用いられる。
【００３７】
このエネルギー関数は、フレームシーケンスの上のアライメント変換をトラッキングするのに役立ち、または、γ→∞の場合、例えば、３つの隣接するフレームを使用して平均的変換を計算するのに役立つ。シングル・フレーム・エネルギーと同様に、変換χfは、変換モデルの複雑さを順次増加させることによって、次第に精製される。
【００３８】
非剛体レジストレーション
グローバル・アフィン変換モデルが、呼吸によって、および、心臓の鼓動によって、生じるインターモダル形状矛盾を完全に補償できるというわけではない。本発明の実施形態は、３Ｄ冠状動脈ツリー中心線と、２つのキャリブレーションされた蛍光透視画像と、の間の視覚的な相関関係を改善できる非剛体レジストレーション方法を提供する。本発明の一実施形態に係る方法は、ユーザが相関関係を識別する必要がないという点で自動といえる。
【００３９】
本発明の一実施形態に係る非剛体変換は、３Ｄ並進ベクトルrs,qのセットrとして表され、３Ｄ並進ベクトルrs,qは、座標（ＣＳ）のＵにおける対応する中心線点xs,pに与えられる。次に、レジストレーション点
【数３９】
は
【数４０】
を用いて計算される。溶液の品質を測定するために規定された、本発明の一実施形態に係るエネルギーは、以下（８）式の形を有する。
【数４１】
本発明の一実施形態に係る最小化プロセスは反復され、位置rs,p=[0,0,0]T|∀{s,q}から開始し、変換rs,qは、勾配降下法に従うことによって徐々に調整される。これは、レジストレーションされた中心線の滑らかな動きを結果として生ずる。図３は、非剛体レジストレーションの進行を示し、図中左から右に、０回、３５０回および１２００回の反復後のサンプル入力画像および中心線位置を示している。
【００４０】
画像エネルギー
画像エネルギーという用語は、再構築された中心線の汎関数である。本発明の一実施形態によると、照合および再構築プロセスを用いて、３Ｄ中心線モデルにおいて全ての点xsqに対する再構築点
【数４２】
を計算することによって、各蛍光透視画像における中心線の射影を構成することができる。図４のフローチャートを参照して、照合は各透視面上で別個に以下のように行われる。本発明の一実施形態に係る照合方法は、ステップ４１において、３Ｄ中心線の制御点を蛍光透視画像nに射影することによって開始する。蛍光透視画像n上の３Ｄ中心線点
【数４３】
の射影x′から開始して、ステップ４２において、蛍光透視画像nに対応するセグメンテーション画像Bn上の検索が、中心線の射影に垂直な方向に行われる。この方向は、各位置において、偏導関数
【数４４】
に基づいて計算される。その結果、ステップ４３において、多くても２つの照合点xa"、xb"は各透視面n∈{1,2}上で識別される。ステップ４１〜４３は、他の蛍光透視画像のためにステップ４４から繰り返される。ステップ４５において、エピポーラ拘束
【数４５】
を最小化するペア(x1,i",x2,i"),i∈a,bが保たれる。ここで、d(x1,i",em,i")は点xn,i"と、xm,i"によって導入されるエピポーラ線との間のユークリッド距離であり、{n,m}∈{1,2}かつn≠mである。ステップ４６において、式（９）を最小化するペアは、システムの射影形状を使用して、再構築３Ｄ点
【数４６】
を計算するために用いられる。ただ１つの点xn,i"が画像面上に見られる場合、それは初期設定により保たれる。図５（ａ）は、中心線x'からの点と２Ｄ構造からの点との照合を示している。点xa"、xb"は、有力候補として識別され、xb"はx'に近いため保たれる。次に、ステップ４１〜４６は、３Ｄ中心線の全制御点のためにステップ４７から繰り返される。
【００４１】
本発明の実施形態では、中心線の射影に局所的に垂直な方向において照合点の検索が行われる。構造のグローバル・アライメントが正しいということは先験的に知られているので、垂直検索方向は、最短点手法と比較して、多数の非現実的な照合を回避する。このことは、２Ｄ自動セグメンテーション（automatic segmentation）が完璧ではない場合に有効である。例えば、３Ｄ射影が、自動的にセグメント化された２Ｄ曲線より長い曲線になった場合、最短点をとることは、チップからの曲線の著しい短縮である。一方、３Ｄ中心線射影のチップにおいて照合は見出されないため、垂直方向での検索は、この効果を低減する。
【００４２】
再構築の品質は保証されず、本発明の実施形態は指標関数
【数４７】
を用いて異常値を棄却する。（１）あらゆる画像面において照合点が見つからない場合、（２）x'とxn"との間の２Ｄ距離がMax2Dより長い場合、（３）xs,qと
【数４８】
との間の距離がMax3Dより長い場合、H(xs,p)=0である。ここで、Max2DおよびMax3Dは、所定の閾値である。実施形態によると、一例として、限定的ではないが、Max2Dは５０ピクセルであり、Max3Dは６ｍｍである。
【００４３】
実施形態によると、式（８）における画像エネルギーは再構築点
【数４９】
とレジストレーションされた中心線点との間の二乗誤差の合計として、全ての有効な点に関して以下のように定義される。
【数５０】
全てのレジストレーションされた中心線点が再構築点と同一の位置にあるとき、このエネルギーは最小となる。
【００４４】
内部エネルギー
本発明の一実施形態によれば、正規化（regularizer）を用いて、視覚的にレジストレーションされた中心線を視覚的にコヒーレントにかつ幾何学的に妥当に保つことができる。本発明の実施形態に従う以下の３つの内部エネルギーの項は、このために考慮可能である。
【数５１】
ここで、μ、νおよびλがエネルギー・バランス・フリー・パラメータである。第１の項EDisp(r)は、非剛体変換による変位が小さいとき最小である。この正規化は、一対の２Ｄ血管造影によって、わずかにしか拘束されない３Ｄ面外動作を最小化する。第２の項ESmooth(r)は、各血管セグメント上の平滑性を確実にすることができる。ここで、
【数５２】
および
【数５３】
は、それぞれ、セグメントパラメータに対する位置(s,q)のための並進ベクトルｒの第１および第２の微分係数を表す。セグメント接合部で、ノイマン境界条件
【数５４】
が仮定される。
【００４５】
はじめの２つの内部エネルギー項については、中心線モデルはフレキシブルであり、接合ノードにおいて以外は、血管の動作は互いに独立している。実際は、血管動作は、血管自身の固有の剛性によってだけでなく、血管が心筋に取り付けられているという理由によっても機械的に拘束される。この点において、第３の項、EMyocard(r)は、心筋制約の最小モデルとしての機能し、血管分岐で剛性の特定の程度を確保し、小さいセグメントがより大きいセグメント上へ倒れこむのを防止できる。本発明の実施形態は、図５（ｂ）に示すように、３つのセグメントの接合部の周りの人工的なリンクを用いて、この制約をモデル化し、図５（ｂ）は、接合部における３つの心筋制約（k1、k2、k3）の形成を示す。Dsを距離パラメータとし、Nに設定される値と接合部に存在する血管の平均直径とを掛け算する。Nの例示的かつ非限定的な値は３である。接合部から距離Ds程離れた各セグメント上のノードは、機械ばねに類似した方法で、他のセグメントの対応するノードに連結され、各リンクは、レジストレーションプロセスの間、拡大または圧縮される場合、その変位比率と比例する力を発生させている。式（１１）において。Kは、ともに結合された全対のノード{i,j}の組として定義される。
【００４６】
エネルギー最小化
本発明の一実施形態による完全なエネルギー関数は、式（８）、式（１０）および式（１１）を使用して定められ、そのオイラー・ラグランジェ方程式を計算して、勾配降下方法に従うことによって、本発明の実施形態において最小化できる。本発明の一実施形態によれば、反復rs,p=[0,0,0]T|∀{s,q}かつ反復z=0での点から開始して、結果として生じる離散更新方程式は、式（１２）のとおりである。
【数５５】
ここで、数ステップεは、小さい定数であり、Ks,q⊂Kは、ノード(s,q)に作用する制約の中のエンプティセットである。本発明の実施形態では、式（１２）は、全体の変位が収束するまで、反復的して評価されうる。より形式的には、全体の変位は、
【数５６】
を意味し、停止基準は、ｚ＞βに関してΔz-Δz-β<ζΔβとして定められる。図６は、本発明の一実施形態に係る、非剛体レジストレーションの間、反復回数の関数として、中心線の全体の変位を示している。曲線の形状から推論されうるように、最終結果はβおよびζの特定の値にあまり影響されない。例示的かつ非限定的な値は、β＝100およびζ＝0.001である。
【００４７】
パラメータ選択
本発明の一実施形態によれば、３つのパラメータμ、υおよびλは、式（１２）の挙動を支配し、特定のタスクに対して選択可能である。特定のパラメータセット（μ＝0.10、υ＝10.0、λ＝1.0）から開始して、以下に詳細に記載されているように、各値を個別に変化させることによって、平均２Ｄ射影エラーおよび３Ｄエラーを計算することによって、方法の挙動についての洞察が得られる。この方法は、後述するように、アフィンアライメントの後、臨床データセットの患者４、および、合成変形を有するデータセットに適用された。前者はＲＣＡを表し、後者はＬＣＡを表す。図８は、非剛体レジストレーション後の平均２Ｄおよび３Ｄエラーを示す。図８において、上の行はパラメータμのデータを示し、中間の行はパラメータυのデータを示し、下の行はパラメータγのデータを示す。各グラフに対して、１つのパラメータが変化し、他の２つのパラメータは一定に保たれた。デフォルト・パラメータ値は、μ＝0.1、υ＝10.0、λ＝1.09である。ｘ軸のスケールは、対数関数である。
【００４８】
パラメータμは全体の変位（displacement）を拘束することができ、２Ｄエラー曲線は計算されたエラーがその値とともに増加することを示している。しかしながら、３Ｄエラー曲線は、μが０．０５から０．１０の範囲の値を有するとき、最小値を示す。加えて、収束率は、μとともに増加する。例えば、レジストレーションを用いて、μ＝{0.001,0.01,0.1}をそれぞれ有する４８８０、３３２０および、１１００の反復に収束した図８（左上）のグラフを生成する。
【００４９】
パラメータυは、モデルの剛性を制御できる。この場合、２Ｄ射影エラーは、０．５０から２．００の範囲の値を使用することを提案する。しかしながら、３Ｄエラー曲線は、より高い値がより正確なレジストレーションになりうることを示す。本発明の実施形態は、例示的かつ非限定的な値として、２Ｄエラー曲線と３Ｄエラー曲線との間の良好なトレードオフを示すυ＝10.0を使用する。
【００５０】
パラメータλは、心筋制約のバランスをとり、グローバル形状の一般的な整合性を保存できる。臨床ＲＣＡデータセットについては、高い値が運動を過度に拘束し、エラーの増加を引き起こすことが判明した。しかしながら、シミュレーションされたＬＣＡデータセットについては、状況はほぼ逆に現れ、１０．０までのλ値は、明らかに有益である。この状況は、ＬＣＡおよびＲＣＡのデータセットに影響を及ぼす変形の性質における違いによって、おそらく説明可能である。例示的かつ非限定的な値λ＝0.5が選択された。なぜなら、それがエラー測定のトレードオフを表し、かつ、図７にて示すように、多くのケースの視覚的外観を改善するためである。図７は、非剛体レジストレーション方法のパラメータλの挙動を、左側のλ＝0および右側のλ＝0.3に関して示している。ＬＣＡおよびＲＣＡデータセットに対して、複数の値のλを用いることができる。本発明の実施形態では、特に言及しない限り、以下の全ての実験において、例示であり限定的ではないパラメータμ=0.5、υ=10.0およびλ=0.5が用いられる。
【００５１】
図９には、本発明の一実施形態に係る非剛体２Ｄ／３Ｄレジストレーション方法のフローチャートを示す。方法は、ステップ９０において、術前の画像体積、例えばＣＴＡ体積からセグメント化された冠動脈ツリーの３Ｄ中心線表現を提供することによって開始する。方法は、ステップ９１において、２つ以上の分割された、冠動脈血管の２Ｄ蛍光透視画像に対する冠動脈ツリーの３Ｄ中心線表現の最初のグローバル・アライメントを実行することによって続行する。実施形態では、上述したように、このグローバル・アライメントは、連続的な並進のみの変換、剛体変換およびアフィン変換によって実行可能である。ステップ９２において、図４を参照して説明したように、３Ｄ中心線の点は、蛍光透視画像における対応する血管中心線点に照合され、照合された２Ｄ点から３Ｄ中心線を再構築する。ステップ９３において、非剛体レジストレーションを実行する。次に、ステップ９３において、エネルギー汎関数を最小化することによって非剛体レジストレーションを実行する。ここで、エネルギー汎関数は、再構築点
【数５７】
とレジストレーションされた中心線点との間の二乗差の総和と、二乗３Ｄレジストレーション・ベクトルの総和と、二乗微分係数３Ｄレジストレーション・ベクトルの総和と、心筋分岐エネルギーと、を含む。
【００５２】
実験結果
本発明の実施形態に係るアルゴリズムは、インターベンション中にレジストレーションされたモデルのライブ蛍光透視画像に対するオーバレイを示すプロトタイピング環境においてＣ＋＋で実施された。Ｃ＋＋での実施は例示であり限定的ではなく、本発明の他の実施形態では、他のコンピュータ言語での実施が可能であることを理解されたい。全てのグローバル・オプティマイザはシングルスレッドである。本発明の一実施形態に係る非剛体レジストレーション法は、OpenMPを使用して部分的にマルチスレッドされている。与えられるランタイムはi7 Q820 Quad Core Intel CPU用である。示されたランタイムは２Ｄ画像セグメンテーションを構成しない。このオペレーションは、画像ごとに平均して０．２００秒のランタイムを有し、標準偏差は０．０８８秒である。シミュレーションと臨床データの両方に関する実験が行われた。
【００５３】
シミュレーション
本発明の実施形態によるアライメントおよびレジストレーションの方法は、シミュレーションのセットにおいてテストされ、いくらか理想的なシナリオにおける本発明の実施形態に従う各種アルゴリズムの性能特性を与え、臨床データを使用する場合には不可能な３Ｄレジストレーションエラーを評価した。中心線は、以下の通り生成され、５１２×５１２ピクセルのサイズおよび０．３４５ピクセル／ｍｍの解像度を有する、一対のデジタル的に再構築されたＸ線画像（ＤＲＲ）にレジストレーションされた。（１）ＣＴボリュームにおいて、セグメントされた冠動脈は、シミュレーションされた透視面に投影され、（２）結果として生じる射影の逆数は、参照背景画像から減算された。図１０（ａ）〜（ｄ）は、ＤＲＲの作成を例示する。ＣＴボリュームにおいて、セグメントされた冠状動脈血栓は、図１０（ａ）のシミュレーションされたＸ線撮影平面に投影される。次に、この画像は、図１０（ｂ）に示すように実際の造影剤のない（contrast-free）血管造影法から減算され、図１０（ｃ）の最終的なＤＲＲを生成する。一実験において、図１０（ｄ）に示すように、σ＝50を有するガウス雑音が加えられる。シミュレーションは、後述する３つのシナリオにおいて使われた。すべての変換がシミュレーション中に公知であるので、本願明細書において提示されるすべての３Ｄエラーは、正確なポイント・ツー・ポイント・エラーである。視覚効果を生じないときでも、エラーのこのスタイルがセグメントの圧縮または拡大を不利にすることに留意する必要がある。それゆえ、それは、大部分のＴＲＥ形成より厳格である。２Ｄエラーは、上述したように、参照セグメンテーションとして実際の中心線の射影を用いて計算される。
【００５４】
（１）初期の解決法に対する依存
中心線は一対のＤＲＲによってレジストレーションされ、異なる初期位置から開始する。これは、初期位置の摂動と、同様に、C、P1、P2またはTの剛体変換行列のミスキャリブレーションに対するグローバル・アライメント方法の感度を定量化する機能を果たす。選択された初期位置は、以下の通りである。T0を装置のキャリブレーションによって、与えられる変換であると仮定する。次に、１２の初期ポイントは３本の主軸に沿って±αｍｍだけ変位したT0と、は３本の主軸の周りに±β度だけ回転したT0である。実施形態は、α∈10×{0,1,...,6}mmおよびβ∈3.75×{0,1,...,6}度を使用する。平方二乗平均エラーは、摂動（並進または回転）の同じレベルを有する全結果に対して計算された。図１１は、初期位置の摂動に関して、オプティマイザの性能を例示するものであり、左列のグラフは剛体アライメントの結果を示し、右列のグラフはアフィンアライメントの結果を示し、上行は方向における変化（変位）の結果を示し、下行は位置における変化の結果を示す。下行において、失った点は、平方二乗平均エラーが３ｍｍ超であることを表す。図１１の結果は、回転摂動が１０度未満のとき、提示される本発明の実施形態に係るすべてのアルゴリズムが満足であることを示す。
【００５５】
（２）画像ノイズに対するロバスト性
分散を有する複数の量のガウス雑音は、蛍光透視画像に加えられ、本発明の実施形態に従う方法の性能上の画質品質の影響を評価した。各オプティマイザは、各ノイズ・レベルに対する１３の初期位置から開始した（上述したように、α＝5mmおよびβ＝7.5度がプラスされたT0）。図１２は、本発明の実施形態に従って、ノイズ・レベルに対する入力ＤＲＲノイズの標準偏差に関して、オプティマイザの性能を例示し、左列のグラフはリジッドアライメントの結果を表し、右列のグラフはアフィンアライメントの結果を表す。σ∈[0,60]を有するガウス白色雑音モデルが用いられた。この試験の結果は、本発明の実施形態による方法が蛍光透視のガウス雑音の存在に影響されないことを示す傾向がある。臨床設定ではほとんど観察されないσ＝60のレベルでさえ、２Ｄおよび３Ｄエラー・レベルは、追加されたノイズの量にわずかな相関しか示さない。これは、選択された２Ｄセグメンテーション・アルゴリズムの性能が良好であるという徴候である。
【００５６】
（３）非剛体変形
これらのシミュレーションにおいて、中心線は、ＣＴで走査した領域をカバーしている３×３×３ノード薄板スプライン（ＴＰＳ）変形モデルを使用して変形した。各ＴＰＳノードは、変形パラメータξに応じた係数によって、中心ノードへ移動する。９つの左のノードは、係数ξによってシフトされ、９つの中心のノードは、0.5によってシフトされ、９つの右のノードは、係数2ξによってシフトされた。図１３は、シミュレーションされたデータを有する若干のサンプルの非剛体レジストレーションを示す。上行は初期位置を示し、下行は最終位置を示す。本発明の一実施形態による非剛体レジストレーション・アルゴリズムは、変形レベルξ∈{0.020,0.025,...,0.060}に対して試験され、中心線の平均３Ｄ変位は、［1.209〜3.638］ｍｍであり、中心線の最大３Ｄ変位は、［2.472〜7.416］ｍｍであった。非変形中心線曲線を、初期の解決法として用い、非剛体レジストレーション・アルゴリズムを適用する前に、グローバル・アライメントは、実行されなかった。図１４は、定量的評価を示しており、シミュレーションされた非剛体変形レベルに関する３Ｄ残差のグラフであり、ベースラインは初期の解決法のエラーに対応する。このアルゴリズムがレベルξ＝0.045までは良好な結果を示すことが判明した。この点では、平均的３Ｄエラーは２．７２１ｍｍから１．１９８ｍｍに減少し、このことは、図１３に示すように質的に重要である。以下、本願明細書に示されるように、性能のこのレベルは臨床シナリオにおいて適切である。
【００５７】
臨床データ
以下のセクションに示される実験に、５つのデータセットが使われた。各データセットは、拡張末期（データセット１、３）または収縮末期（データセット２、４、５）にて得られる１つのＣＴＡ走査と、１つの２方向Ｘ線透視法記録と、を含む。２つのモダリティはＥＣＧゲーティングを用いて時間的にアライメントされ、冠動脈は専門家によって、ＣＴＡにおいて、半自動式にセグメントされた。セグメントされた中心線の平均３Ｄインターポイント距離は、［1.36〜1.60］ｍｍの範囲にあった。血管造影図は、３つの場合の左冠状動脈（ＬＣＡ）と、その他２の場合の右冠動脈（ＲＣＡ）と、を撮像する。画像サイズは全ての場合において、５１２×５１２であり、画像解像度はデータセット１〜５に対してそれぞれ｛0.345,0.279,0.279,0.216,0.279｝ｍｍ／ピクセルである。獲得率は、１５ｆｐｓである。標準Ｃアームキャリブレーション行列が使われた。必須ではないが、キャリブレーション行列は、データセットごとに異なるにもかかわらず、獲得の間一定であった。ＣＴＡおよび蛍光透視試験は、冠動脈疾患の治療のための患者に定められた。
【００５８】
本発明の一実施形態によるアライメントおよびレジストレーション・アルゴリズムの性能の評価は、平均２Ｄ射影エラー
【数５８】
を使用して定量化された。ここで、Mは３Ｄ中心線モデルの点の総数であり、
【数５９】
は参照２ＤセグメンテーションSnの距離変換であり、対応する蛍光透視画像上で手動追跡によって得られた。
【００５９】
グローバル・アライメント：オプティマイザの性能評価
本発明の実施形態は、グローバル・アライメントのために、７つのローカル（Ｌ）アルゴリズムおよび２つのグローバル（Ｇ）アルゴリズムを検査した。以下、これらのアルゴリズムを簡潔に説明する。
【００６０】
ベスト・ネイバー（Best Neighbor）（Ｌ）
各反復で、パラメータ空間のすべての主要な方向におけるステップが考慮される。最善の移動が適用され、新規の反復が開始する。移動が解決法を改善しないとき、ステップは半分にされる。
【００６１】
ネルダー・ミード（別名滑降シンプレックス）（Ｌ）
一般のシンプレックス（またはポリエイラクブカ）のｍ＋１頂点で関数値を評価することによって、ｍ次元関数を最小化する古典的な数の最適化方法である。各反復で、最悪の値を有する頂点は、拡大動作または縮小動作の一方が続く反転動作を使用して、別のものに置換される。ネルダー等による「A Simplex Method For Function Minimization」、Ｃｏｍｐｕｔ．Ｊ．第７巻、第４号、３０８〜３１３頁、１９６５年参照。
【００６２】
Ｓｂｐｌｘ３（Ｌ）
この方法は、低次元の部分空間の問題を分解し、ネルダー・ミード・アルゴリズムを使用して、検索を実行する。オンラインで利用可能なＳ．Ｇ．ジョンソンによる「The Nlopt Nonlinear-Optimization Package」（http://ab-initio.mit.edu/nlopt）参照。
【００６３】
Ｃｏｂｙｌａ（Ｌ）
線形近似による位相制約下最適化（Constrained Optimization By Linear Approximation）方法は、頂点シンプレックスを用いたコスト関数および制約の線形近似を構成することによって、機能し、この近似値を最小化する。シンプレックスの半径は次第に減少し、標準形状を維持する。パウエルによる「Advances in Optimization and Numerical Analysis」（A Direct Search Optimization Method That Models the Objective and Constraint Functions by Linear Interpolation、ドルドレヒト、クルーワー アカデミック、１９９４年、５１―６７頁）、および、パウエルによる「Direct Search Algorithms For Optimization Calculations」（Dept. Appl. Math. Theoretical Phys., Cambridge Univ., U.K., Tech. Rep., 1998）参照。
【００６４】
Ｂｏｂｙｑａ（Ｌ）
二字近似による境界最適化（Bound Optimization By Quadratic Approximation）方法の各反復は、概して展補間点を考慮することによって、構成されたコスト関数の二次近似を使用する。パウエルによる「The BOBYQA Algorithm For Bound Constrained Optimization Without Derivatives」、（Centre Math. Sci., Cambridge Univ., Cambridge, U.K., Tech. Rep., 2009）参照。
【００６５】
パウエル・ブレント（Ｌ）
各反復で、一連の正確な１Ｄ線最適化は、ブレントの方法を使用して実行される。その方法は、共役（conjugate）検索方向のパウエルの方法を使用して更新される。ブレントによる「Minimization Without Derivatives, Englewood Cliffs, NJ: Prentice-Hall, 1973」、および、パウエルによる「An Efficient Method For Finding The Minimum Of A Function Of Several Variables Without Calculating Derivatives」、Comput. J.、第７巻、第２号、１５５〜１６２頁、１９６４年参照。
【００６６】
Ｐｒａｘｉｓ（Ｌ）
ブレントのＰＲｉｎｃｉｐａｌ―ＡＸＩＳ方法は、共役検索方向のパウエルの方法を更新したものであり、「Algorithms for Minimization Without Derivatives, Englewood Cliffs, NJ: Prentice-Hall, 1973」に記載されている。
【００６７】
微分展開（Ｇ）
これは、集団ベースの確率的な大域的最適化方法であり、その主要な特徴は、集団進化の基礎としての個人の対の違いのベクトルの使用である。ストーン等による「Differential Evolution-A Simple And Efficient Heuristic For Global Optimization Over Continuous Spaces」、J. Global Optimizat.、第１１巻、３４１〜３５９頁、１９９７年参照。
【００６８】
Ｄｉｒｅｃｔ（Ｇ）
この大域的最適化アルゴリズムは、本明細書の場合のように、有限境界拘束に関する課題のために設計される。パラメータ空間は、ますます小さいハイパーレクタングルに分割することによって、系統的かつ決定論的に検索される。ジョーンズ等による「Lipschitzian Optimization Without The Lipschitz Constant」、Optimizat. Theory Appl.、第７９巻、１５７〜１８１頁、１９９３年参照。
【００６９】
各アルゴリズムの自由パラメータは、アルゴリズムの原作者またはメーカーの勧告に従って調整された。以下のアルゴリズムの実施は、ＳＢＰＬＸ、Ｃｏｂｙｌａ、ＢｏｂｙｑａおよびＤｉｒｅｃｔについて、http://ab-initio.mit.edu/nloptのＮＬｏｐｔライブラリから得られる。課題の性質を変えないために、多数の境界は、オプティマイザに供給され、トランスレーションパラメータ用の±２００ｍｍ、回転パラメータ用の±４５度、スケール・パラメータ用の［0.5,1.5］、スケール・回転パラメータ用の±３６０度である[式（２）においてA1012]。
【００７０】
微分展開を除いて、基準を停止する同一セットが、全オプティマイザ用に用いられた。具体的には、オプティマイザは、１つの最適化の後、（１）コスト関数の値が、1e-14未満で減少するか、あるいは、（２）最適化されたパラメータの値が、並進パラメータに関して1e-4ｍｍ未満、回転パラメータに関して1e-5ラジアン、スケーリング・パラメータに関して1e-7だけ変更するまで、進むことが許可された。パラメータ許容度は、直観的な解釈を有する。例えば、1e-4ｍｍの並進の違いは、画像平面上の１／２０００ピクセルの最大変位に対応する。微分展開オプティマイザに関しては、確率的なアルゴリズムの収束が規則的ではないので、本発明の実施形態は、コスト関数を、約１０秒の中央値の動作に対応する１０００００回評価するまで、オプティマイザが動作するのを許容した。
【００７１】
式（１）にて表される本発明の実施形態に従うグローバル・アライメント・エネルギーを最小化する際および良質の２Ｄ／３Ｄアライメントを生成する際の９つのオプティマイザの性能は、以下の実験的な設定を使用して評価された。５人の患者のデータセットの各々に対して、３対の蛍光透視画像が考慮された。それは、ＣＴＡ獲得に時間的にアライメントされた画像および前後の隣接するフレームである。時間的アライメントは２つのモダリティ間の観察可能な違いを最小化するのを支援し、それによって、レジストレーションに失敗する危険を減らす。レジストレーション行列が毎秒ほぼ１回更新され、外科的ガイダンスを改善するため、これは臨床的に合理的である。加えて、オプティマイザのロバスト性を試験するために、１３の異なる初期点が、初期化のために使われた。T0を、装置のキャリブレーションによって与えられる変換であると仮定する。１３の初期点は、T0と、３本の主軸に沿って±５ｍｍだけ位置がずれたT0と、３本の主軸の周りに±７．５度だけ回転したT0と、である。合計５×３×１３＝１９５回の実験は、テストされる各オプティマイザに対して実行された。上述したように、式（１）は、（１）並進のみの変換、（２）剛体変換および（３）アフィン変換を連続して用いることにより、最小化可能である。
【００７２】
全９つのオプティマイザに対する３つの変換モデルの各々によって、最小化した後残った平均残余のエネルギーの平方二乗平均の分析は、アフィン変換モデルを有するエネルギー関数を最小化するための最善のアルゴリズムが以下の、微分展開、ネルダー・ミード、パウエル・ブレント、ベスト・ネイバー、ＳｂｐｌｘおよびＤｉｒｅｃｔであるということを明らかにする。平均２Ｄ射影エラーが考慮されると、順番はわずかに変化して、微分展開、パウエル・ブレント、ベスト・ネイバー、Ｓｂｐｌｘ、ネルダー・ミードおよびＤｉｒｅｃｔとなる。違いは、アライメントプロセスで使用される自動セグメンテーションと、エラーを計算するために使用される手動セグメンテーションと、の間の相違におそらく起因している。それにもかかわらず、上部にランクを付けされたアルゴリズムの性能は、満足であるように見える。ＢｏｂｙｑａおよびＣｏｂｙｌａアルゴリズムの比較的低い性能は、それぞれ、エネルギー関数の形状がそれらの二次モデルおよび線形モデルによって、良好に表されないという徴候であるかもしれない。Ｐｒａｘｉｓ方法は、この課題に最も不適切なローカル・オプティマイザであるように見える。
【００７３】
微分展開アルゴリズムが、残余エネルギーおよび平均２Ｄエラーの観点で最高の性能を有すると判明した。しかしながら、平均２Ｄエラーの増加は、ローカル・アルゴリズムより長い約１００回の計算時間に重要ではなかった。しかし、そのオプティマイザの計算時間の中央値は、ちょうど１１秒未満であり、タイム・クリティカルでないアプリケーションには適当となりうる。一般に、グローバル・オプティマイザは、最も高性能のローカル・オプティマイザより低いエネルギーまたはエラー形状に、極めてまれに至るということに留意されたい。これは、後者がグローバルな最適条件の近くにある解決法を実際に見つけるという徴候である。グローバル・オプティマイザが最低の解決法を回避することに一般に成功する点に注意する価値がある。残余エネルギーと認められた視覚の対応との間に良好な相関が存在する。それにもかかわらず、アフィン変換モデルで最高のアライメントに関してさえ、比較的大きい相違が、投影中心線と２Ｄ蛍光透視との間に存在する。これは、モダリティ間の時間的アライメントの不正確性による、データセットの非アフィン変形の存在と、さまざまな取得ステップの間、患者の位置の変化により導入される形状変化と、におそらく起因している。
【００７４】
グローバル・アライメント方法と非剛体レジストレーションとの比較
ベスト・ネイバー・オプティマイザを用いて得られた詳細なグローバル・アライメントの結果が表示され、非剛体レジストレーション後に得られた結果と比較される。全ての場合において、中心線モデルは、ＣＴ再構築用に使用される同一の心臓位相において、ゲートされる蛍光透視フレームでレジストレーションされた。並進、剛体、アフィンおよび非剛体レジストレーション後の結果は、図３および図１５の３人の患者からのデータセットに関して示され、図１５では、左から右に、入力画像、剛体アライメント、アフィンアライメントおよび非剛体レジストレーションを示す。両方の２方向画像は、すべての場合に示される。矢印は関心領域を示す。２Ｄマニュアル・セグメンテーションに関する２つの透視面上の中心線の平均残余の２Ｄ射影エラーも、５人の患者に対する剛体、アフィンおよび非剛体レジストレーション後に、（ミリメータで）算出され、その結果は図１６の表Ｉに示されている。図１５の上の行および下の行は、それぞれ、表Ｉの患者１および患者２に対応する。図１５の２つの左下の画像フレームのＣＴＯラベルは、慢性完全閉塞を意味する。非剛体レジストレーションに必要な計算時間は、７３０ｍｓ〜３３００ｍｓの範囲にあり、平均１５９１ｍｓであった。
【００７５】
モデルの複雑さが増加するにつれて、平均２Ｄ射影エラーおよびそれらの標準偏差は減少するか、または、ほぼ一定のままである。しかしながら、剛体、アフィンおよび非剛体変換モデルの相対的な貢献は、データセットごとに変化する。これは、個々のケースにより示される変形の性質によって説明可能であり、変形の性質は、患者の位置、ＣＴＡと２Ｄ蛍光透視獲得との間の間隔および獲得プロトコルに関連している。また、患者の呼吸によって、非剛体の心臓の変形が生じうる。心臓の鼓動は非剛体変形の原因でもあるため、ＥＣＧゲーティングによる時間的アライメントの精度も重要である。最大エラーもモデル複雑さに伴って減少する傾向があるが、時々安定なままである。このような場合、自動のセグメンテーション方法により計算されるため、レジストレーション・アルゴリズムの性能は不完全な２Ｄセグメンテーションにより制限されるかもしれない。実際、ほとんど情報が利用できない領域では、本発明の実施形態による非剛体レジストレーション方法は、中心線をできるだけ変えてはならない。この特性は、慢性完全閉塞（ＣＴＯ）の場合のように、構造の一部がＸ線蛍光透視画像上で十分には見えない状況において望ましい。例えば、図１５の下行は、中心線がＣＴＯの上下に良好にアライメントされているが、コントラストが低い部分にほとんど変形がないことを示す。
【００７６】
マルチ・フレーム・シナリオのグローバル・アライメント
２方向の一対に対する中心線のアライメントの他の潜在的に役立つシナリオは、２方向血管造影法からの時間的フレームシーケンスに対する中心線のアライメントである。マルチ・フレーム・シナリオにおける本発明の一実施形態によるグローバル・アライメント方法の性能は、患者１のＬＣＡデータセットから一連の１１の２方向フレームを使用して示される。グローバル・アライメントは、本発明の実施形態に従う剛体およびアフィン変換モデルを使用して計算された。式（５）のフレーム間組織化パラメータγは、０．０〜１．０の間で変化した。式（５）のコスト関数を最小化するのに最も効果的であると判明したので、ベスト・ネイバー・オプティマイザがこの場合使われた。得られたサンプル結果は図１７に示され、図１７はフレームシーケンスの上のレジストレーションを示す。図１７において、一番上の行は、γ＝0.1およびEMulti*＝120.284を有する剛体変換モデルを表し、中間の行は、γ＝0およびEMulti*＝119.356を有するアフィン変換モデルを表し、一番下の行は、γ＝0.1およびEMulti*＝116.905を有する正規化アフィン変換モデルを表し、ここで、EMulti*は、式（５）の右辺の第１項に対応する。矢印は関心領域を示す。計算時間は、一番上の行、中間の行および一番下の行で示される実験に対して、それぞれ、３９秒、８９秒および９１秒であった。予想されるように、アフィン変換モデルを使用すると剛体変換モデルを使用するより良好な定性的結果に至ることが判明した。アフィン変換モデルを使用するとき、組織化（γ＝０．１）を使用することも結果を改善した。大きい組織化値（γ＞１．０）を用いると、モデルの剛性は増加し、オプティマイザは初期位置から遠い極小値を見つけることができない。予備調査は、例えば、パラメータが協調して変化可能なことによって、領域に特有のオプティマイザがより良好な結果を達成しうることを示唆する。それにもかかわらず、本発明の実施形態に従う設定は、大部分の心周期の間、ＬＣＡのアライメント変換をトラッキングするのに適切に聞こえた。ＲＣＡデータセットに関する実験は、説得力のある結果につながらなかった。事実、本発明の一実施形態によるグローバル・アライメント手順が機能しても、ＣＴＡおよびＸ線取得が近い時間的アライメントに近い場合を除き、心周期の間、ＲＣＡにより支持される剛体変形の量によって、アフィン・モデルの効果がなくなる。
【００７７】
右冠動脈の半自動のトラッキング
この実験では、本発明の実施形態による非剛体レジストレーション方法に基づく半自動手順を用いて、ＲＣＡのトラッキングを行った。ゲートされたフレームから開始して、本発明の一実施形態によるグローバル・アライメントおよび非剛体レジストレーションが実行される。次に、本発明の一実施形態によるこの変形した中心線モデルは、次の一対の蛍光透視フレームのための第一のモデルおよび位置としてオペレータにより用いられる。図１８に示されるように、このプロセスは１つの心周期上の全フレームのために繰り返され、図１８は、左から右に、ゲートされたフレーム、＋４、＋６、＋８および＋１０フレームを表す。第一の曲線は１８１であり、レジストレーションされた曲線は１８２である。了解されるように、この手順は、１つの心周期の間、ＲＣＡのトラッキングを成功させる。この実験がシミュレーション研究の前に行われたので、わずかに異なるパラメータセット、μ＝０．０５、υ＝５．０およびλ＝０．１が用いられた。必要な時間は約１５分であり、これはオフラインの設定において、合理的である。
【００７８】
考察
臨床データの実験に基づいて、ネルダー・ミード、パウエル・ブレント、ベスト・ネイバーおよびＳｂｐｌｘのローカル・オプティマイザが短時間（中央値＜１０５ｍｓ、最大値＜９５０ｍｓ）で良好な結果を与える。シミュレーションも考慮すると、ネルダー・ミード・アルゴリズムは、全体的に最善の性能を与えた。２つのグローバル・オプティマイザは、その結果が大きな改良につながることは極めてまれであり、一桁高い計算時間を必要とした。計算時間を無視すると、微分展開アルゴリズムは、一般に、最善の解決案を提供した。ローカル・オプティマイザを使用するとき、アライメント時間は一貫して１秒未満であるため、この方法は、インターベンション中の使用に適するものになる。剛体変換の代わりにアフィン変換を使用することの利点は、データセットの性質に依存し、場合によっては重要になりうる。３Ｄおよび２ＤモダリティはＥＣＧゲーティングを使用して時間的にアライメントされるとき、全体として、グローバル・アライメント手順がＬＣＡおよびＲＣＡデータセットの両方の使用に適切であること判明した。ＣＴＯ患者からのデータセットに関する実験は、類似の臨床ケースに適用されるとき、本発明の実施形態に従う方法の利点を示す。本発明の実施形態による方法は、ＣＴＡ獲得からセグメント化された３Ｄ形状を有する２Ｄ画像を増加することによって、インターベンション・ガイダンスを支援することもでき、このことにより２Ｄ画像の解釈に固有の曖昧性を著しく低減する。
【００７９】
非剛体レジストレーションに関しては、本発明の一実施例による変形モデルに含まれる正規化は、中心線を最小量だけ変形させることによって、残りの構造との一貫性を確実に保つことによって、欠測データ（例えば、図１５の患者２）を有する領域の管理を可能にする。過剰な変形または矯正は発生しない。全体の計算時間は、通常、３秒未満であり、これはインタラクティブなアプリケーションに使用できる。本発明の一実施形態による非剛体レジストレーション方法は、適切な初期点から開始するとき、隣接する２Ｄ構造に対する中心線をスナップし、特に困難なローコントラストＣＴＯケースに関して、２Ｄ画像の解釈をより容易にする。これも、操作上の画像と一緒に更新された３Ｄモデルを示すことを可能にし、３Ｄ空間の改良された認識を提供できる。
【００８０】
最終的に、複数のＸ線２方向血管造影法フレームを使用する実験によって、本発明の一実施形態によるマルチフレーム・グローバル・アライメント方法が、ＬＣＡデータセットに良好に作用するとともに、本発明の一実施形態によるフレーム間の正規化が、改良された結果につながるということが判明した。本発明の一実施形態による非剛体半自動のトラッキング方法は、より大きな非剛体変形を有するケース、例えば非時間的にアライメントされたＲＣＡデータセットを有する場合を扱い、そのようなデータセットに適用された。本発明の一実施形態による半自動の方法は、１つの心拍をカバーしているシーケンスに対して約１５分の相互作用を必要とし、それは実用的なオフラインのアプリケーションに合理的に思える。
【００８１】
システムの実現
本発明の実施形態は、各種形態のハードウェア、ソフトウェア、ファームウェア、特殊目的のプロセスあるいはその組み合わせによって実現可能であることを理解されたい。一実施形態では、本発明は、コンピュータ可読のプログラム記憶装置に具体的に組み込まれたアプリケーションプログラムとしてソフトウェアで実現可能である。アプリケーションプログラムは、任意の適切な構造を有する機械にアップロード可能であるとともに、当該機械によって実行可能である。
【００８２】
図１９は、ライブ蛍光透視画像を用いた冠動脈モデルの非剛体２Ｄ／３Ｄレジストレーションを実現するコンピュータシステムのブロック図である。図１９を参照すると、本発明を実施するためのコンピュータシステム191は、とりわけ、中央処理装置（ＣＰＵ）192と、メモリ193と、入出力（Ｉ／Ｏ）インタフェース194と、を具えることができる。コンピュータシステム191は、一般的に、Ｉ／Ｏインタフェース194を介してディスプレイ195と、マウスやキーボードのような複数の入力装置196と、に接続されている。支持回路は、キャッシュ、電源、クロック回路および通信バスのような回路を含むことができる。メモリ193は、ランダムアクセスメモリ（ＲＡＭ）、読み出し専用メモリ（ＲＯＭ）、ディスクドライブ、テープドライブ等、あるいはその組み合わせを含むことができる。本発明は、ルーティン197として実現され、このルーティン197はメモリ193に記憶され、ＣＰＵ192によって実行され、信号源198からの信号を処理する。それゆえ、コンピュータシステム191は、汎用コンピュータシステムであり、本発明のルーティン197を実行するとき特殊目的のコンピュータシステムになる。
【００８３】
コンピュータシステム191は、また、オペレーティングシステムと、マイクロ命令コードと、を含む。本明細書に記載した様々な処理および機能は、オペレーティングシステムを介して実行されるマイクロ命令コードの一部あるいはアプリケーションプログラムの一部（あるいはその組み合わせ）とすることができる。さらに、様々なその他の周辺装置、例えば、追加のデータ記憶装置および印刷装置を、コンピュータプラットフォームに接続することができる。
【００８４】
添付図面に示されたシステムの構成要素および方法のステップのいくつかは、ソフトウェアによって実施可能であるため、システムの構成要素（あるいはプロセスステップ）間の実際の接続は、本発明がプログラムされる方法に依存して異なるということを理解されたい。本明細書に記載した本発明の教示を考慮すると、当業者は本発明と同一および類似の実施形態あるいは構成を考えることができる。
【００８５】
本発明は、例示的実施形態を参照して詳細に説明してきたが、当業者であれば、添付の特許請求の範囲に規定した本発明の精神および範囲から逸脱することなく各種変更および置換が可能であることを理解するであろう。
【特許請求の範囲】
【請求項１】
心臓インターベンション中の２Ｄ蛍光透視画像を用いたデジタル３Ｄ冠動脈モデルの非剛体レジストレーションのための方法であって、前記方法は、
Qsの３Ｄ制御点からなる一組のSセグメントを具える３Ｄ中心線であって、冠動脈ツリーのデジタル化された３Ｄ中心線の表示を提供するステップと、
前記３Ｄ中心線を、少なくとも２つの２Ｄ蛍光透視画像に、グローバル・アライメントするステップと、
エネルギー汎関数を最小化することによって、前記３Ｄ中心線を、少なくとも２つの２Ｄ蛍光透視画像に、非剛体レジストレーションするステップと、
を含み、
前記エネルギー汎関数は、再構築された中心線点とレジストレーションされた中心線点との間の二乗差の総和と、二乗３Ｄレジストレーション・ベクトルの総和と、二乗微分係数３Ｄレジストレーション・ベクトルの総和と、心筋分岐エネルギーと、を含み、
前記３Ｄ中心線の非剛体レジストレーションは、前記３Ｄ中心線の座標系にて、対応する中心線点xs,qに適用される一組の３Ｄ並進ベクトルrs,qとして表される、
ことを特徴とする方法。
【請求項２】
前記冠動脈ツリーは前記２Ｄ蛍光透視画像からセグメントされ、前記セグメントされた冠動脈ツリーは前記２Ｄ蛍光透視画像の各々の２値化として表される、
請求項１に記載の方法。
【請求項３】
前記３Ｄ中心線を、前記少なくとも２つの２Ｄ蛍光透視画像に、グローバル・アライメントするステップは、
前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線の並進のみのレジストレーションを実行するステップと、
前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線の剛体レジストレーションを実行するステップと、
前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線のアフィン・レジストレーションを実行するステップと、
を含み、
各ステップでは、
【数１】
の形のエネルギー汎関数は最小化され、
ｊは、前記並進のみのレジストレーションと、前記剛体レジストレーションと、前記アフィン・レジストレーションと、を分類し、
Cは、前記３Ｄ中心線の座標系と３Ｄ画像形成装置の座標系との間の座標変換であり、
Tjは、前記並進のみのレジストレーションと、前記剛体レジストレーションと、前記アフィン・レジストレーションと、に対して、前記３Ｄ画像形成装置の座標系と、前記蛍光透視画像を取得する２方向Ｃアームの座標系と、の間の座標変換であり、
xは、均一な座標の３Ｄ点であり、
nは、前記蛍光透視画像上の合計を表し、
sおよびqは、前記３Ｄ中心線の前記セグメント上の制御点を分類し、
Pnは、前記２方向Ｃアームの前記座標系と、第ｎ番目の蛍光透視画像の前記座標系と、の間の座標変換であり、
Ψnは、前記３Ｄ位置ｘを前記第ｎ番目の蛍光透視画像に射影する射影オペレータであり、
【数２】
は、前記第ｎ番目の蛍光透視画像の２Ｄ点と、前記第ｎ番目の二値化された画像Bnの前記セグメントされた冠動脈ツリーに最も近い点と、の間の距離である、
請求項２に記載の方法。
【請求項４】
変換行列Tjは、再センタリングされ、
前記変換行列Tjの回転パラメータの変化は、前記３Ｄ画像形成装置の前記座標系の原点とは対照的に、前記３Ｄ中心線の前記座標系の原点の周りで発生するように見える運動を誘導する、
請求項３に記載の方法。
【請求項５】
時間的に連続したF対の２Ｄ蛍光透視画像を提供するステップと、
前記エネルギー汎関数
【数３】
を最小化することによって、前記３Ｄ中心線を、各対の２Ｄ蛍光透視画像にグローバル・アライメントするステップと、
をさらに含み、
fは時点の合計であり、
Tfは各時点の変換であり、
γは所定のパラメータであり、
D(Tf・C・xs,q,Tf+1・C・xs,q)はフレーム間距離である、
請求項３に記載の方法。
【請求項６】
前記３Ｄ中心線を少なくとも２つの２Ｄ蛍光透視画像にグローバル・アライメントするステップの後、
前記３Ｄ中心線の点を、前記２Ｄ蛍光透視画像の対応する血管中心線点に照合し、前記照合した２Ｄ点から３Ｄ中心線点
【数４】
を再構築するステップをさらに含む、
請求項１に記載の方法。
【請求項７】
前記点照合は、
各３Ｄ中心線点xs,pの射影x'を、前記グローバル・アライメントからの各２Ｄ蛍光透視画像上に見つけるステップと、
x'から開始して、各２Ｄ蛍光透視画像に対応するセグメンテーション画像Bnを、前記セグメンテーション画像B上の前記３Ｄ中心線の射影に垂直な線に沿って検索し、中心線点xs,pに照合する２つまでの点を見つけるステップと、
xs,pに最も近い点を照合点xnとして選択するステップと、
前記蛍光透視画像の照合点を用いて、前記心臓インターベンションに用いられている画像処理システムの射影形状を使用している３Ｄ中心線点
【数５】
を再構築するステップと、
を含む、
請求項６に記載の方法。
【請求項８】
再構築された中心線点
【数６】
とレジストレーションされた中心線点
【数７】
との二乗差の前記総和は、
【数８】
であり、
【数９】
は指標関数
【数１０】
であり、
照合点がいかなる画像平面でも見つからない場合、あるいは、２Ｄ蛍光透視画像上の
【数１１】
の射影x'と前記第ｎ番目の２Ｄ蛍光透視画像の照合点xn"との間の２Ｄ距離が所定の閾値Max2Dより大きい場合、あるいは、中心線点xs,qと対応する再構築された点
【数１２】
との間の３Ｄ距離が所定の閾値Max3Dより大きい場合、H(xs,p)=0であり、
その他の場合、H(xs,p)=1である、
請求項７に記載の方法。
【請求項９】
二乗３Ｄレジストレーション・ベクトルの前記総和と、二乗微分係数３Ｄレジストレーション・ベクトルの前記総和と、前記エネルギー汎関数の前記心筋分岐エネルギーは、
【数１３】
であり、
rは、レジストレーション・ベクトルであり、
xは、前記３Ｄ中心線上の制御点であり、
sおよびqの前記総和は、前記３Ｄ中心線の全セグメントの全制御点の合計であり、
Kは、前記３Ｄ中心線の接合部から距離Ds離れた各セグメントの点を、前記接合部で接合される他のセグメントの対応する制御点に連結することにより形成される制御点の全対の一組であり、
前記接合部では、連結がレジストレーションの間拡大または圧縮される場合、各連結はその変位に比例した力を発生させる、
請求項８に記載の方法。
【請求項１０】
前記エネルギー汎関数は、
【数１４】
を反復することによって、最小化され、
zは、点rs,p=[0,0,0]T|∀{s,q}から反復z=0で開始する反復カウンタであり、
εは、小さい定数であり、
Ks,q⊂Kは、制御点(s,q)に作用する一組の制約であり、
μ、υおよびλは、所定のパラメータである、
請求項９に記載の方法。
【請求項１１】
心臓インターベンション中の２Ｄ蛍光透視画像を用いたデジタル３Ｄ冠動脈モデルの非剛体レジストレーションのための方法であって、前記方法は、
Qsの３Ｄ制御点からなる一組のSセグメントを具える３Ｄ中心線であって、冠動脈ツリーのデジタル化された３Ｄ中心線の表示を提供するステップと、
前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線の並進のみのレジストレーションを実行することによって、前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線の剛体レジストレーションを実行することによって、前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線のアフィン・レジストレーションを実行することによって、前記３Ｄ中心線を少なくとも２つの２Ｄ蛍光透視画像にグローバル・アライメントするステップと、
を含み、
各ステップでは、
【数１５】
の形のエネルギー汎関数は最小化され、
ｊは、前記並進のみのレジストレーションと、前記剛体レジストレーションと、前記アフィン・レジストレーションと、を分類し、
Cは、前記３Ｄ中心線の座標系と３Ｄ画像形成装置の座標系との間の座標変換であり、
Tjは、前記並進のみのレジストレーションと、前記剛体レジストレーションと、前記アフィン・レジストレーションと、に対して、前記３Ｄ画像形成装置の座標系と、前記蛍光透視画像を取得する２方向Ｃアームの座標系と、の間の座標変換であり、
xは、均一な座標の３Ｄ点であり、
nは、前記蛍光透視画像上の合計を表し、
sおよびqは、前記３Ｄ中心線の前記セグメント上の制御点を分類し、
Pnは、前記２方向Ｃアームの前記座標系と、第ｎ番目の蛍光透視画像の前記座標系と、の間の座標変換であり、
Ψnは、前記３Ｄ位置xを前記第ｎ番目の蛍光透視画像に射影する射影オペレータであり、
【数１６】
は、前記第ｎ番目の蛍光透視画像の２Ｄ点と、前記第ｎ番目の二値化された画像Bnの前記セグメントされた冠動脈ツリーに最も近い点と、の間の距離であり、
前記方法は、前記３Ｄ中心線を、前記少なくとも２つの２Ｄ蛍光透視画像に、非剛体レジストレーションするステップを含む、
ことを特徴とする方法。
【請求項１２】
前記３Ｄ中心線を、前記少なくとも２つの２Ｄ蛍光透視画像に、非剛体レジストレーションするステップは、
前記３Ｄ中心線の点を、前記２Ｄ蛍光透視画像の対応する血管中心線点に照合し、前記照合した２Ｄ点から３Ｄ中心線点
【数１７】
を再構築するステップと、
エネルギー汎関数を最小化するステップと、
を含み、
前記エネルギー汎関数は、再構築された中心線点とレジストレーションされた中心線点との間の二乗差の総和と、二乗３Ｄレジストレーション・ベクトルの総和と、二乗微分係数３Ｄレジストレーション・ベクトルの総和と、心筋分岐エネルギーと、を含み、
前記３Ｄ中心線の非剛体レジストレーションは、前記３Ｄ中心線の座標系にて、対応する中心線点xs,qに適用される一組の３Ｄ並進ベクトルrs,qとして表される、
請求項１１に記載の方法。
【請求項１３】
前記エネルギー汎関数は、
【数１８】
であり、
rは、レジストレーション・ベクトルであり、
xは、前記３Ｄ中心線上の制御点であり、
【数１９】
は、再構築された中心線点であり、
【数２０】
は、レジストレーションされた中心線点であり、
【数２１】
は指標関数
【数２２】
であり、
照合点がいかなる画像平面でも見つからない場合、あるいは、２Ｄ蛍光透視画像上の
【数２３】
の射影x'と前記第ｎ番目の２Ｄ蛍光透視画像の照合点xn"との間の２Ｄ距離が所定の閾値Max2Dより大きい場合、あるいは、中心線点xs,qと対応する再構築された点
【数２４】
との間の３Ｄ距離が所定の閾値Max3Dより大きい場合、H(xs,p)=0であり、
その他の場合、H(xs,p)=1であり、
sおよびqの前記総和は、前記３Ｄ中心線の全セグメントの全制御点の合計であり、
Kは、前記３Ｄ中心線の接合部から距離Ds離れた各セグメントの点を、前記接合部で接合される他のセグメントの対応する制御点に連結することにより形成される制御点の全対の一組である、
請求項１２に記載の方法。
【請求項１４】
コンピュータ可読の非一時的プログラム記憶装置であって、
前記コンピュータにより実行されるとともに、心臓インターベンション中の２Ｄ蛍光透視画像を用いたデジタル３Ｄ冠動脈モデルの非剛体レジストレーションのための方法を実行する命令プログラムを有し、前記方法は、
Qsの３Ｄ制御点からなる一組のSセグメントを具える３Ｄ中心線であって、冠動脈ツリーのデジタル化された３Ｄ中心線の表示を提供するステップと、
前記３Ｄ中心線を、少なくとも２つの２Ｄ蛍光透視画像に、グローバル・アライメントするステップと、
エネルギー汎関数を最小化することによって、前記３Ｄ中心線を、少なくとも２つの２Ｄ蛍光透視画像に、非剛体レジストレーションするステップと、
を含み、
前記エネルギー汎関数は、再構築された中心線点とレジストレーションされた中心線点との間の二乗差の総和と、二乗３Ｄレジストレーション・ベクトルの総和と、二乗微分係数３Ｄレジストレーション・ベクトルの総和と、心筋分岐エネルギーと、を含み、
前記３Ｄ中心線の非剛体レジストレーションは、前記３Ｄ中心線の座標系にて、対応する中心線点xs,qに適用される一組の３Ｄ並進ベクトルrs,qとして表される、
ことを特徴とするプログラム記憶装置。
【請求項１５】
前記冠動脈ツリーは前記２Ｄ蛍光透視画像からセグメントされ、前記セグメントされた冠動脈ツリーは前記２Ｄ蛍光透視画像の各々の２値化として表される、
請求項１４に記載のコンピュータ可読のプログラム記憶装置。
【請求項１６】
前記３Ｄ中心線を、前記少なくとも２つの２Ｄ蛍光透視画像に、グローバル・アライメントするステップは、
前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線の並進のみのレジストレーションを実行するステップと、
前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線の剛体レジストレーションを実行するステップと、
前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線のアフィン・レジストレーションを実行するステップと、
を含み、
各ステップでは、
【数２５】
の形のエネルギー汎関数は最小化され、
ｊは、前記並進のみのレジストレーションと、前記剛体レジストレーションと、前記アフィン・レジストレーションと、を分類し、
Cは、前記３Ｄ中心線の座標系と３Ｄ画像形成装置の座標系との間の座標変換であり、
Tjは、前記並進のみのレジストレーションと、前記剛体レジストレーションと、前記アフィン・レジストレーションと、に対して、前記３Ｄ画像形成装置の座標系と、前記蛍光透視画像を取得する２方向Ｃアームの座標系と、の間の座標変換であり、
xは、均一な座標の３Ｄ点であり、
nは、前記蛍光透視画像上の合計を表し、
sおよびqは、前記３Ｄ中心線の前記セグメント上の制御点を分類し、
Pnは、前記２方向Ｃアームの前記座標系と、第ｎ番目の蛍光透視画像の前記座標系と、の間の座標変換であり、
Ψnは、前記３Ｄ位置ｘを前記第ｎ番目の蛍光透視画像に射影する射影オペレータであり、
【数２６】
は、前記第ｎ番目の蛍光透視画像の２Ｄ点と、前記第ｎ番目の二値化された画像Bnの前記セグメントされた冠動脈ツリーに最も近い点と、の間の距離である、
請求項１５に記載のコンピュータ可読のプログラム記憶装置。
【請求項１７】
変換行列Tjは、再センタリングされ、
前記変換行列Tjの回転パラメータの変化は、前記３Ｄ画像形成装置の前記座標系の原点とは対照的に、前記３Ｄ中心線の前記座標系の原点の周りで発生するように見える運動を誘導する、
請求項１６に記載のコンピュータ可読のプログラム記憶装置。
【請求項１８】
前記方法は、
時間的に連続したF対の２Ｄ蛍光透視画像を提供するステップと、
前記エネルギー汎関数
【数２７】
を最小化することによって、前記３Ｄ中心線を、各対の２Ｄ蛍光透視画像にグローバル・アライメントするステップと、
をさらに含み、
fは、時点の合計であり、
Tfは、各時点の変換であり、
γは所定のパラメータであり、
D(Tf・C・xs,q,Tf+1・C・xs,q)は、フレーム間距離である、
請求項１６に記載のコンピュータ可読のプログラム記憶装置。
【請求項１９】
前記方法は、
前記３Ｄ中心線を少なくとも２つの２Ｄ蛍光透視画像にグローバル・アライメントするステップの後、
前記３Ｄ中心線の点を、前記２Ｄ蛍光透視画像の対応する血管中心線点に照合し、前記照合した２Ｄ点から３Ｄ中心線点
【数２８】
を再構築するステップをさらに含む、
請求項１４に記載のコンピュータ可読のプログラム記憶装置。
【請求項２０】
前記点照合は、
各３Ｄ中心線点xs,pの射影x'を、前記グローバル・アライメントからの各２Ｄ蛍光透視画像上に見つけるステップと、
x'から開始して、各２Ｄ蛍光透視画像に対応するセグメンテーション画像Bnを、前記セグメンテーション画像B上の前記３Ｄ中心線の射影に垂直な線に沿って検索し、中心線点xs,pに照合する２つまでの点を見つけるステップと、
xs,pに最も近い点を照合点xnとして選択するステップと、
前記蛍光透視画像の照合点を用いて、前記心臓インターベンションに用いられている画像処理システムの射影形状を使用している３Ｄ中心線点
【数２９】
を再構築するステップと、
を含む、
請求項１９に記載のコンピュータ可読のプログラム記憶装置。
【請求項２１】
再構築された中心線点
【数３０】
とレジストレーションされた中心線点
【数３１】
との二乗差の総和は、
【数３２】
であり、
【数３３】
は指標関数
【数３４】
であり、
照合点がいかなる画像平面でも見つからない場合、あるいは、２Ｄ蛍光透視画像上の
【数３５】
の射影x'と前記第ｎ番目の２Ｄ蛍光透視画像の照合点xn"との間の２Ｄ距離が所定の閾値Max2Dより大きい場合、あるいは、中心線点xs,qと対応する再構築された点
【数３６】
との間の３Ｄ距離が所定の閾値Max3Dより大きい場合、H(xs,p)=0であり、
その他の場合、H(xs,p)=1である、
請求項２０に記載のコンピュータ可読のプログラム記憶装置。
【請求項２２】
二乗３Ｄレジストレーション・ベクトルの前記総和と、二乗微分係数３Ｄレジストレーション・ベクトルの前記総和と、前記エネルギー汎関数の前記心筋分岐エネルギーは、
【数３７】
であり、
rは、レジストレーション・ベクトルであり、
xは、前記３Ｄ中心線上の制御点であり、
sおよびqの前記総和は、前記３Ｄ中心線の全セグメントの全制御点の合計であり、
Kは、前記３Ｄ中心線の接合部から距離Ds離れた各セグメントの点を、前記接合部で接合される他のセグメントの対応する制御点に連結することにより形成される制御点の全対の一組であり、
前記接合部では、連結がレジストレーションの間拡大または圧縮される場合、各連結はその変位に比例した力を発生させる、
請求項２１に記載のコンピュータ可読のプログラム記憶装置。
【請求項２３】
前記エネルギー汎関数は、
【数３８】
を反復することによって、最小化され、
zは、点rs,p=[0,0,0]T|∀{s,q}から反復z=0で開始する反復カウンタであり、
εは、小さい定数であり、
Ks,q⊂Kは、制御点(s,q)に作用する一組の制約であり、
μ、υおよびλは、所定のパラメータである、
請求項２２に記載のコンピュータ可読のプログラム記憶装置。
【請求項１】
心臓インターベンション中の２Ｄ蛍光透視画像を用いたデジタル３Ｄ冠動脈モデルの非剛体レジストレーションのための方法であって、前記方法は、
Qsの３Ｄ制御点からなる一組のSセグメントを具える３Ｄ中心線であって、冠動脈ツリーのデジタル化された３Ｄ中心線の表示を提供するステップと、
前記３Ｄ中心線を、少なくとも２つの２Ｄ蛍光透視画像に、グローバル・アライメントするステップと、
エネルギー汎関数を最小化することによって、前記３Ｄ中心線を、少なくとも２つの２Ｄ蛍光透視画像に、非剛体レジストレーションするステップと、
を含み、
前記エネルギー汎関数は、再構築された中心線点とレジストレーションされた中心線点との間の二乗差の総和と、二乗３Ｄレジストレーション・ベクトルの総和と、二乗微分係数３Ｄレジストレーション・ベクトルの総和と、心筋分岐エネルギーと、を含み、
前記３Ｄ中心線の非剛体レジストレーションは、前記３Ｄ中心線の座標系にて、対応する中心線点xs,qに適用される一組の３Ｄ並進ベクトルrs,qとして表される、
ことを特徴とする方法。
【請求項２】
前記冠動脈ツリーは前記２Ｄ蛍光透視画像からセグメントされ、前記セグメントされた冠動脈ツリーは前記２Ｄ蛍光透視画像の各々の２値化として表される、
請求項１に記載の方法。
【請求項３】
前記３Ｄ中心線を、前記少なくとも２つの２Ｄ蛍光透視画像に、グローバル・アライメントするステップは、
前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線の並進のみのレジストレーションを実行するステップと、
前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線の剛体レジストレーションを実行するステップと、
前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線のアフィン・レジストレーションを実行するステップと、
を含み、
各ステップでは、
【数１】
の形のエネルギー汎関数は最小化され、
ｊは、前記並進のみのレジストレーションと、前記剛体レジストレーションと、前記アフィン・レジストレーションと、を分類し、
Cは、前記３Ｄ中心線の座標系と３Ｄ画像形成装置の座標系との間の座標変換であり、
Tjは、前記並進のみのレジストレーションと、前記剛体レジストレーションと、前記アフィン・レジストレーションと、に対して、前記３Ｄ画像形成装置の座標系と、前記蛍光透視画像を取得する２方向Ｃアームの座標系と、の間の座標変換であり、
xは、均一な座標の３Ｄ点であり、
nは、前記蛍光透視画像上の合計を表し、
sおよびqは、前記３Ｄ中心線の前記セグメント上の制御点を分類し、
Pnは、前記２方向Ｃアームの前記座標系と、第ｎ番目の蛍光透視画像の前記座標系と、の間の座標変換であり、
Ψnは、前記３Ｄ位置ｘを前記第ｎ番目の蛍光透視画像に射影する射影オペレータであり、
【数２】
は、前記第ｎ番目の蛍光透視画像の２Ｄ点と、前記第ｎ番目の二値化された画像Bnの前記セグメントされた冠動脈ツリーに最も近い点と、の間の距離である、
請求項２に記載の方法。
【請求項４】
変換行列Tjは、再センタリングされ、
前記変換行列Tjの回転パラメータの変化は、前記３Ｄ画像形成装置の前記座標系の原点とは対照的に、前記３Ｄ中心線の前記座標系の原点の周りで発生するように見える運動を誘導する、
請求項３に記載の方法。
【請求項５】
時間的に連続したF対の２Ｄ蛍光透視画像を提供するステップと、
前記エネルギー汎関数
【数３】
を最小化することによって、前記３Ｄ中心線を、各対の２Ｄ蛍光透視画像にグローバル・アライメントするステップと、
をさらに含み、
fは時点の合計であり、
Tfは各時点の変換であり、
γは所定のパラメータであり、
D(Tf・C・xs,q,Tf+1・C・xs,q)はフレーム間距離である、
請求項３に記載の方法。
【請求項６】
前記３Ｄ中心線を少なくとも２つの２Ｄ蛍光透視画像にグローバル・アライメントするステップの後、
前記３Ｄ中心線の点を、前記２Ｄ蛍光透視画像の対応する血管中心線点に照合し、前記照合した２Ｄ点から３Ｄ中心線点
【数４】
を再構築するステップをさらに含む、
請求項１に記載の方法。
【請求項７】
前記点照合は、
各３Ｄ中心線点xs,pの射影x'を、前記グローバル・アライメントからの各２Ｄ蛍光透視画像上に見つけるステップと、
x'から開始して、各２Ｄ蛍光透視画像に対応するセグメンテーション画像Bnを、前記セグメンテーション画像B上の前記３Ｄ中心線の射影に垂直な線に沿って検索し、中心線点xs,pに照合する２つまでの点を見つけるステップと、
xs,pに最も近い点を照合点xnとして選択するステップと、
前記蛍光透視画像の照合点を用いて、前記心臓インターベンションに用いられている画像処理システムの射影形状を使用している３Ｄ中心線点
【数５】
を再構築するステップと、
を含む、
請求項６に記載の方法。
【請求項８】
再構築された中心線点
【数６】
とレジストレーションされた中心線点
【数７】
との二乗差の前記総和は、
【数８】
であり、
【数９】
は指標関数
【数１０】
であり、
照合点がいかなる画像平面でも見つからない場合、あるいは、２Ｄ蛍光透視画像上の
【数１１】
の射影x'と前記第ｎ番目の２Ｄ蛍光透視画像の照合点xn"との間の２Ｄ距離が所定の閾値Max2Dより大きい場合、あるいは、中心線点xs,qと対応する再構築された点
【数１２】
との間の３Ｄ距離が所定の閾値Max3Dより大きい場合、H(xs,p)=0であり、
その他の場合、H(xs,p)=1である、
請求項７に記載の方法。
【請求項９】
二乗３Ｄレジストレーション・ベクトルの前記総和と、二乗微分係数３Ｄレジストレーション・ベクトルの前記総和と、前記エネルギー汎関数の前記心筋分岐エネルギーは、
【数１３】
であり、
rは、レジストレーション・ベクトルであり、
xは、前記３Ｄ中心線上の制御点であり、
sおよびqの前記総和は、前記３Ｄ中心線の全セグメントの全制御点の合計であり、
Kは、前記３Ｄ中心線の接合部から距離Ds離れた各セグメントの点を、前記接合部で接合される他のセグメントの対応する制御点に連結することにより形成される制御点の全対の一組であり、
前記接合部では、連結がレジストレーションの間拡大または圧縮される場合、各連結はその変位に比例した力を発生させる、
請求項８に記載の方法。
【請求項１０】
前記エネルギー汎関数は、
【数１４】
を反復することによって、最小化され、
zは、点rs,p=[0,0,0]T|∀{s,q}から反復z=0で開始する反復カウンタであり、
εは、小さい定数であり、
Ks,q⊂Kは、制御点(s,q)に作用する一組の制約であり、
μ、υおよびλは、所定のパラメータである、
請求項９に記載の方法。
【請求項１１】
心臓インターベンション中の２Ｄ蛍光透視画像を用いたデジタル３Ｄ冠動脈モデルの非剛体レジストレーションのための方法であって、前記方法は、
Qsの３Ｄ制御点からなる一組のSセグメントを具える３Ｄ中心線であって、冠動脈ツリーのデジタル化された３Ｄ中心線の表示を提供するステップと、
前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線の並進のみのレジストレーションを実行することによって、前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線の剛体レジストレーションを実行することによって、前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線のアフィン・レジストレーションを実行することによって、前記３Ｄ中心線を少なくとも２つの２Ｄ蛍光透視画像にグローバル・アライメントするステップと、
を含み、
各ステップでは、
【数１５】
の形のエネルギー汎関数は最小化され、
ｊは、前記並進のみのレジストレーションと、前記剛体レジストレーションと、前記アフィン・レジストレーションと、を分類し、
Cは、前記３Ｄ中心線の座標系と３Ｄ画像形成装置の座標系との間の座標変換であり、
Tjは、前記並進のみのレジストレーションと、前記剛体レジストレーションと、前記アフィン・レジストレーションと、に対して、前記３Ｄ画像形成装置の座標系と、前記蛍光透視画像を取得する２方向Ｃアームの座標系と、の間の座標変換であり、
xは、均一な座標の３Ｄ点であり、
nは、前記蛍光透視画像上の合計を表し、
sおよびqは、前記３Ｄ中心線の前記セグメント上の制御点を分類し、
Pnは、前記２方向Ｃアームの前記座標系と、第ｎ番目の蛍光透視画像の前記座標系と、の間の座標変換であり、
Ψnは、前記３Ｄ位置xを前記第ｎ番目の蛍光透視画像に射影する射影オペレータであり、
【数１６】
は、前記第ｎ番目の蛍光透視画像の２Ｄ点と、前記第ｎ番目の二値化された画像Bnの前記セグメントされた冠動脈ツリーに最も近い点と、の間の距離であり、
前記方法は、前記３Ｄ中心線を、前記少なくとも２つの２Ｄ蛍光透視画像に、非剛体レジストレーションするステップを含む、
ことを特徴とする方法。
【請求項１２】
前記３Ｄ中心線を、前記少なくとも２つの２Ｄ蛍光透視画像に、非剛体レジストレーションするステップは、
前記３Ｄ中心線の点を、前記２Ｄ蛍光透視画像の対応する血管中心線点に照合し、前記照合した２Ｄ点から３Ｄ中心線点
【数１７】
を再構築するステップと、
エネルギー汎関数を最小化するステップと、
を含み、
前記エネルギー汎関数は、再構築された中心線点とレジストレーションされた中心線点との間の二乗差の総和と、二乗３Ｄレジストレーション・ベクトルの総和と、二乗微分係数３Ｄレジストレーション・ベクトルの総和と、心筋分岐エネルギーと、を含み、
前記３Ｄ中心線の非剛体レジストレーションは、前記３Ｄ中心線の座標系にて、対応する中心線点xs,qに適用される一組の３Ｄ並進ベクトルrs,qとして表される、
請求項１１に記載の方法。
【請求項１３】
前記エネルギー汎関数は、
【数１８】
であり、
rは、レジストレーション・ベクトルであり、
xは、前記３Ｄ中心線上の制御点であり、
【数１９】
は、再構築された中心線点であり、
【数２０】
は、レジストレーションされた中心線点であり、
【数２１】
は指標関数
【数２２】
であり、
照合点がいかなる画像平面でも見つからない場合、あるいは、２Ｄ蛍光透視画像上の
【数２３】
の射影x'と前記第ｎ番目の２Ｄ蛍光透視画像の照合点xn"との間の２Ｄ距離が所定の閾値Max2Dより大きい場合、あるいは、中心線点xs,qと対応する再構築された点
【数２４】
との間の３Ｄ距離が所定の閾値Max3Dより大きい場合、H(xs,p)=0であり、
その他の場合、H(xs,p)=1であり、
sおよびqの前記総和は、前記３Ｄ中心線の全セグメントの全制御点の合計であり、
Kは、前記３Ｄ中心線の接合部から距離Ds離れた各セグメントの点を、前記接合部で接合される他のセグメントの対応する制御点に連結することにより形成される制御点の全対の一組である、
請求項１２に記載の方法。
【請求項１４】
コンピュータ可読の非一時的プログラム記憶装置であって、
前記コンピュータにより実行されるとともに、心臓インターベンション中の２Ｄ蛍光透視画像を用いたデジタル３Ｄ冠動脈モデルの非剛体レジストレーションのための方法を実行する命令プログラムを有し、前記方法は、
Qsの３Ｄ制御点からなる一組のSセグメントを具える３Ｄ中心線であって、冠動脈ツリーのデジタル化された３Ｄ中心線の表示を提供するステップと、
前記３Ｄ中心線を、少なくとも２つの２Ｄ蛍光透視画像に、グローバル・アライメントするステップと、
エネルギー汎関数を最小化することによって、前記３Ｄ中心線を、少なくとも２つの２Ｄ蛍光透視画像に、非剛体レジストレーションするステップと、
を含み、
前記エネルギー汎関数は、再構築された中心線点とレジストレーションされた中心線点との間の二乗差の総和と、二乗３Ｄレジストレーション・ベクトルの総和と、二乗微分係数３Ｄレジストレーション・ベクトルの総和と、心筋分岐エネルギーと、を含み、
前記３Ｄ中心線の非剛体レジストレーションは、前記３Ｄ中心線の座標系にて、対応する中心線点xs,qに適用される一組の３Ｄ並進ベクトルrs,qとして表される、
ことを特徴とするプログラム記憶装置。
【請求項１５】
前記冠動脈ツリーは前記２Ｄ蛍光透視画像からセグメントされ、前記セグメントされた冠動脈ツリーは前記２Ｄ蛍光透視画像の各々の２値化として表される、
請求項１４に記載のコンピュータ可読のプログラム記憶装置。
【請求項１６】
前記３Ｄ中心線を、前記少なくとも２つの２Ｄ蛍光透視画像に、グローバル・アライメントするステップは、
前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線の並進のみのレジストレーションを実行するステップと、
前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線の剛体レジストレーションを実行するステップと、
前記少なくとも２つの２Ｄ蛍光透視画像に対する前記３Ｄ中心線のアフィン・レジストレーションを実行するステップと、
を含み、
各ステップでは、
【数２５】
の形のエネルギー汎関数は最小化され、
ｊは、前記並進のみのレジストレーションと、前記剛体レジストレーションと、前記アフィン・レジストレーションと、を分類し、
Cは、前記３Ｄ中心線の座標系と３Ｄ画像形成装置の座標系との間の座標変換であり、
Tjは、前記並進のみのレジストレーションと、前記剛体レジストレーションと、前記アフィン・レジストレーションと、に対して、前記３Ｄ画像形成装置の座標系と、前記蛍光透視画像を取得する２方向Ｃアームの座標系と、の間の座標変換であり、
xは、均一な座標の３Ｄ点であり、
nは、前記蛍光透視画像上の合計を表し、
sおよびqは、前記３Ｄ中心線の前記セグメント上の制御点を分類し、
Pnは、前記２方向Ｃアームの前記座標系と、第ｎ番目の蛍光透視画像の前記座標系と、の間の座標変換であり、
Ψnは、前記３Ｄ位置ｘを前記第ｎ番目の蛍光透視画像に射影する射影オペレータであり、
【数２６】
は、前記第ｎ番目の蛍光透視画像の２Ｄ点と、前記第ｎ番目の二値化された画像Bnの前記セグメントされた冠動脈ツリーに最も近い点と、の間の距離である、
請求項１５に記載のコンピュータ可読のプログラム記憶装置。
【請求項１７】
変換行列Tjは、再センタリングされ、
前記変換行列Tjの回転パラメータの変化は、前記３Ｄ画像形成装置の前記座標系の原点とは対照的に、前記３Ｄ中心線の前記座標系の原点の周りで発生するように見える運動を誘導する、
請求項１６に記載のコンピュータ可読のプログラム記憶装置。
【請求項１８】
前記方法は、
時間的に連続したF対の２Ｄ蛍光透視画像を提供するステップと、
前記エネルギー汎関数
【数２７】
を最小化することによって、前記３Ｄ中心線を、各対の２Ｄ蛍光透視画像にグローバル・アライメントするステップと、
をさらに含み、
fは、時点の合計であり、
Tfは、各時点の変換であり、
γは所定のパラメータであり、
D(Tf・C・xs,q,Tf+1・C・xs,q)は、フレーム間距離である、
請求項１６に記載のコンピュータ可読のプログラム記憶装置。
【請求項１９】
前記方法は、
前記３Ｄ中心線を少なくとも２つの２Ｄ蛍光透視画像にグローバル・アライメントするステップの後、
前記３Ｄ中心線の点を、前記２Ｄ蛍光透視画像の対応する血管中心線点に照合し、前記照合した２Ｄ点から３Ｄ中心線点
【数２８】
を再構築するステップをさらに含む、
請求項１４に記載のコンピュータ可読のプログラム記憶装置。
【請求項２０】
前記点照合は、
各３Ｄ中心線点xs,pの射影x'を、前記グローバル・アライメントからの各２Ｄ蛍光透視画像上に見つけるステップと、
x'から開始して、各２Ｄ蛍光透視画像に対応するセグメンテーション画像Bnを、前記セグメンテーション画像B上の前記３Ｄ中心線の射影に垂直な線に沿って検索し、中心線点xs,pに照合する２つまでの点を見つけるステップと、
xs,pに最も近い点を照合点xnとして選択するステップと、
前記蛍光透視画像の照合点を用いて、前記心臓インターベンションに用いられている画像処理システムの射影形状を使用している３Ｄ中心線点
【数２９】
を再構築するステップと、
を含む、
請求項１９に記載のコンピュータ可読のプログラム記憶装置。
【請求項２１】
再構築された中心線点
【数３０】
とレジストレーションされた中心線点
【数３１】
との二乗差の総和は、
【数３２】
であり、
【数３３】
は指標関数
【数３４】
であり、
照合点がいかなる画像平面でも見つからない場合、あるいは、２Ｄ蛍光透視画像上の
【数３５】
の射影x'と前記第ｎ番目の２Ｄ蛍光透視画像の照合点xn"との間の２Ｄ距離が所定の閾値Max2Dより大きい場合、あるいは、中心線点xs,qと対応する再構築された点
【数３６】
との間の３Ｄ距離が所定の閾値Max3Dより大きい場合、H(xs,p)=0であり、
その他の場合、H(xs,p)=1である、
請求項２０に記載のコンピュータ可読のプログラム記憶装置。
【請求項２２】
二乗３Ｄレジストレーション・ベクトルの前記総和と、二乗微分係数３Ｄレジストレーション・ベクトルの前記総和と、前記エネルギー汎関数の前記心筋分岐エネルギーは、
【数３７】
であり、
rは、レジストレーション・ベクトルであり、
xは、前記３Ｄ中心線上の制御点であり、
sおよびqの前記総和は、前記３Ｄ中心線の全セグメントの全制御点の合計であり、
Kは、前記３Ｄ中心線の接合部から距離Ds離れた各セグメントの点を、前記接合部で接合される他のセグメントの対応する制御点に連結することにより形成される制御点の全対の一組であり、
前記接合部では、連結がレジストレーションの間拡大または圧縮される場合、各連結はその変位に比例した力を発生させる、
請求項２１に記載のコンピュータ可読のプログラム記憶装置。
【請求項２３】
前記エネルギー汎関数は、
【数３８】
を反復することによって、最小化され、
zは、点rs,p=[0,0,0]T|∀{s,q}から反復z=0で開始する反復カウンタであり、
εは、小さい定数であり、
Ks,q⊂Kは、制御点(s,q)に作用する一組の制約であり、
μ、υおよびλは、所定のパラメータである、
請求項２２に記載のコンピュータ可読のプログラム記憶装置。
【図１】
【図４】
【図５】
【図９】
【図１６】
【図１９】
【図２】
【図３】
【図６】
【図７】
【図８】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１７】
【図１８】
【図４】
【図５】
【図９】
【図１６】
【図１９】
【図２】
【図３】
【図６】
【図７】
【図８】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１７】
【図１８】
【公開番号】特開２０１３−７１０１６（Ｐ２０１３−７１０１６Ａ）
【公開日】平成２５年４月２２日（２０１３．４．２２）
【国際特許分類】
【外国語出願】
【出願番号】特願２０１２−２１７５３７（Ｐ２０１２−２１７５３７）
【出願日】平成２４年９月２８日（２０１２．９．２８）
【出願人】（５１０２９９４４４）シーメンス コーポレーション (3)
【氏名又は名称原語表記】Ｓｉｅｍｅｎｓ Ｃｏｒｐｏｒａｔｉｏｎ
【住所又は居所原語表記】１７０ Ｗｏｏｄ Ａｖｅｎｕｅ Ｓｏｕｔｈ， Ｉｓｅｌｉｎ， Ｎｅｗ Ｊｅｒｓｅｙ ０８８３０， Ｕｎｉｔｅｄ Ｓｔａｔｅｓ ｏｆ Ａｍｅｒｉｃａ
【Ｆターム（参考）】
【公開日】平成２５年４月２２日（２０１３．４．２２）
【国際特許分類】
【出願番号】特願２０１２−２１７５３７（Ｐ２０１２−２１７５３７）
【出願日】平成２４年９月２８日（２０１２．９．２８）
【出願人】（５１０２９９４４４）シーメンス コーポレーション (3)
【氏名又は名称原語表記】Ｓｉｅｍｅｎｓ Ｃｏｒｐｏｒａｔｉｏｎ
【住所又は居所原語表記】１７０ Ｗｏｏｄ Ａｖｅｎｕｅ Ｓｏｕｔｈ， Ｉｓｅｌｉｎ， Ｎｅｗ Ｊｅｒｓｅｙ ０８８３０， Ｕｎｉｔｅｄ Ｓｔａｔｅｓ ｏｆ Ａｍｅｒｉｃａ
【Ｆターム（参考）】
[ Back to top ]