先行画像で制約する連続画像再構成の方法
下取り採取した画像データから高品質画像を再構成する方法を提供する。画像再構成方法は、幾つかの異なる画像様式に適用可能である。特に、本発明は、適切な先行画像を逐次画像再構成処理に組込み、前のタイムフレームからの画像フレームを利用し、現在の画像フレームの再構成を抑制する画像再構成方法を提供する。
【発明の詳細な説明】
【技術分野】
【０００１】
本願は、表題"Method for Image Reconstruction Using Prior Image Constrained Compressed Sensing"で出願(2008年1月14日付)している米国仮出願(番号61/020,847)、及び、表題"Method for Image Reconstruction Using Prior Image Constrained Compressed Sensing"で(2008年6月9日付)出願している米国仮出願(番号61/059,891)において優先権を主張している。
【０００２】
この発明は、米国政府機関の国立衛生研究所の支援(NIH EB005712及びNIH EB007021)でなされた。米国政府は、本発明に特定の権利を有する。
【背景技術】
【０００３】
発明の領域は、医用画像に関し、特に画像を再構成する方法に関する。
【０００４】
コンピュータ断層撮影装置において、Ｘ線源は、直交座標系のＸ−Ｙ面（画像面）内で平行になる円錐形状ビームを放射する。Ｘ線ビームは、撮像する被写体（例、内科患者）を通過し、配列した放射線検出器に突き当たる。透過放射線強度は、被写体のＸ線ビームの減衰に依存し、各検出器は、ビーム減衰の測定となる個別の電気信号を生ずる。全ての検出器からの減衰測定は、個別に取得され、透過像、減衰像、又は、投影と呼ばれるものを形成する。
【０００５】
従来のＣＴ装置の光源及び検出器配列は、架台上で画像面内及び被写体周囲で回転させ、 Ｘ線ビームの被写体に交差する角度を常に変化させる。任意の角度での検出器配列からの透過像は、１つの視野となり、被写体の走査は、Ｘ線源及び検出器の一回転での異なる角度方向で形成される一組の視野を含む。２Ｄ走査の場合、データは、被写体から得られる２次元薄片に対応する画像の構成に処理される。２Ｄデータから画像を再構成する一般的な方法は、フィルタ補正逆投影法と称されている。この画像再構成処理は、走査中に取得される減衰測定を、表示部の対応するピクセルの輝度制御に使用されるＣＴ数又はＣＴ値（Hounsfield units）と呼ばれる整数値に変換される。
【発明の概要】
【発明が解決しようとする課題】
【０００６】
本発明は、画像再構成の方法を提供するものである。
【課題を解決するための手段】
【０００７】
本発明の画像再構成法は、正確な再構成の利点を信号対ノイズ比（ＳＮＲ）のゲインと組合せている。つまり、被検者の適切な先行画像を画像再構成処理に組み込んでいる。ＳＮＲの増加に加えて、本発明の方法は、分解能を一次的に増す。更に、得られる所望画像及び投影の数に対し、本発明の方法は、従来の画像再構成法（例、圧縮センシング：ＣＳ）が失敗している正確な画像再構成を与える。
【０００８】
本発明は、時間分解画像データから良質の画像フレームを実時間で再構成する方法を提供する。この方法は、実時間画像方法を利用する画像誘導治療処置に適用できる。例えば、Ｘ線蛍光誘導治療処置は、本発明の方法で利点を得る。この場合、良質の画像フレームが、患者への全Ｘ線量の減少で再構成ができる。
【０００９】
本発明の方法は、時系列の画像フレームでの特定画像フレームの再構成を与える。本方法は、時系列の最も近い画像フレームから現在の画像フレームを差し引く。例えば、前画像フレームが、現在の画像フレームから差し引かれる。このように、各新規に再構成された画像フレームが、動的な先行画像として使用され、現在の画像フレームを疎にする。現在の画像フレームは、更に少ない投影データから正確な再構成を可能にする。
【００１０】
本発明の利点を下記明細書で説明する。記載において、添付図を用いて発明の実施形態を説明している。実施形態は、発明の範囲を限定するものでなく、発明の解釈の参照とするものである。
【図面の簡単な説明】
【００１１】
【図１】本発明で利用している画像再構成法の実施形態の工程図である。
【図２】本発明で利用している画像再構成法の他の実施形態の工程図である。
【図３】本発明で利用している画像再構成法の別の実施形態の工程図である。
【図４】（Ａ）本発明を使用するＣアームＸ線装置の図。（Ｂ）図４ＡのＣアームＸ線装置の構成図である。
【図５】（Ａ）図４ＡのＣアームＸ線装置のＸ線源及び検出器の図。（Ｂ）図４ＡのＣアームＸ線装置で用いられるＣアーム走査経路の図である。
【図６】図６は、本発明を使用する磁気共鳴画像（ＭＲＩ）装置の構成図である。
【図７】図７は、本発明を使用する他の磁気共鳴画像（ＭＲＩ）装置の構成図である。
【図８】図８は、本発明の実施形態の工程を説明する工程図である。
【発明を実施するための形態】
【００１２】
一般的に、一組のデータから画像を再構成する方法は、測定データサンプルＹから所望の画像Ｉを推測する一連の計算ステップを含む。更に具体的には、画像再構成は、次の適合性条件を満たす。
【００１３】
式（１） ＡＩ＝Ｙ
【００１４】
Ａはシステム行列で、一般的には、所望画像Ｉを取得データサンプルＹに関連付ける前方投影演算子として見なせる。コンピュータ断層撮影（ＣＴ）画像の処理の場合、システム行列は再投影操作を含む。式（１）の適合性条件は、別の言葉で言えば、画像が忠実に再構成される時、前方演算は、測定した投影データの正確な推定値を生成するために、実データ取得手続に略似ているべきである。
【００１５】
本発明の実施で使用する画像再構成法は、良質の画像を提供する。「先行画像」は、圧縮センシング（ＣＳ）原理が利用されている反復画像再構成法を抑制するために使用されている。例えば、ＣＳで一般に使用されているスパース変換に加えて、画像を所望画像から先行画像を取り去って更に疎にする。その結果、ＣＳ法で要求されるよりも少ないサンプル数の使用で、画像が正確に再構成される。
【００１６】
更に具体的には、先行画像Ｉｐ及び再構成する所望画像Ｉにおいて、画像再構成の本発明の方法は、次の目的関数の最小化で達成される。
【００１７】
式（２） α||Ψ１(Ｉ−Ｉｐ)||１＋(１−α)||Ψ２Ｉ||1
Ψ１、Ψ２：スパース変換
||...||１：Ｌ１ノルム演算
α：目的関数式（２）の２つの項の相対的重みの制御に利用される正則化変数
【００１８】
上記から、次式
式（３） ||x||１＝Σｉ＝１Ｎ|xｉ|
は、Ｎ次元ベクトルｘのＬ１ノルムを意味する。一般的には、真のＬ１ノルムからの偏差は、所望画像に適切な画質を維持する限り可能である。例えば、目的関数式（２）は、次式に一般化できる。
【００１９】
式（４） α||Ψ１(Ｉ−Ｉｐ)||ｐ＋(１−α)||Ψ２Ｉ||ｐ
||...||ｐ：次の形式を有するＬｐノルム演算
式（５） ||ｘ||ｐ＝(Σｉ＝１Ｎ|ｘｊ|ｐ)１/ｐ
【００２０】
上記のように、好ましくはｐ＝１．０であるが、しかし、異なるｐ値が可能である。
【００２１】
当業者には周知であるが、ｐ値が１．０から更に外れると、再構成した所望画像の劣化が一層明白になる。
【００２２】
式（２）でのスパース変換Ψ１、Ψ２は一般的には異なるが、同一であっても良い。例次的なスパース変換は、ウェーブレット変換、１次差分、２次差分、及び、次式の離散階調変換を含んでいる。
式（６） ▽m,nＩ(m,n)＝√([I(m+1,n)−I(m,n)]２＋[I(m,n＋1)−Ｉ(m,n)]２)
指数ｍ、ｎ：画像Ｉのピクセルの位置
▽m,nＩ(m,n)で特定される画像は、通常、「階調画像」と称されている。
【００２３】
目的関数式（２）の両方の項が重要である。重要性の結果として、正則化変数の選択αが、全体の画像再構成処理の制御に利用される。正則化変数αの選択は、先行画像Ｉｐの選択及び手元の医療用途の選択にも依存する。例えば、目的関数式（２）の第２項は、先行画像Ｉｐから潜在的に受け継ぐ縞のノイズを軽減する。実例として、α≒０．３−０．７の正則化変数の選択は、呼吸器官の画像適用には一般的に十分である。
【００２４】
適合性条件式（１）を全体の画像再構成に一層良く取り込むために、ラグランジュ乗数法が利用される。この場合、適合性条件は、式（２）で設定した目的関数の最小化に更に制約を加えて使用される。新しい目的関数は、次式の形式となる。
式（７） α||Ψ１(Ｉ−Ｉｐ)||1＋(１−α)||Ψ２Ｉ||1＋λ||Ｘ||２２
λ：ラグランジュ乗数
Ｘ：差分行列
||...||２２：平方Ｌ２ノルム演算
Ｎ次元ベクトルｘに対し、Ｌ２ノルムは次式となる。
式（８） ||ｘ||２２＝Σｉ＝１Ｎｘｉ２
【００２５】
式（７）の差分行列は、適合性条件式（１）から成り、次式となる。
式（９） Ｘ＝ＡＩ−Ｙ
【００２６】
ラグランジュ乗数λは、本発明の実施で使用される特定の画像装置に対し経験的に決定される。例えば、ラグランジュ乗数λは、所望のデータ適合性条件と先行画像Ｉｐの類似性との間で所定の交換条件で決定される。大きいラグランジュ乗数λが選択されると、再構成された画像は、低ノイズ分散を有する。しかし、先行画像に特徴な高い空間分解能の損失となる。小さいラグランジュ乗数λが使用されると、先行画像の特徴である高い空間分解能が維持されるが、ノイズ分散が所望の画像で高くなる。このような状況は、利用される画像装置で達成されるコントラスト対ノイズ比に影響を及ぼす。
【００２７】
式（７）に示す目的関数は、画像装置のノイズを説明するために変更できる。この場合、次の目的関数が最小化される。
【００２８】
式（１０） α||Ψ１(Ｉ−Ｉｐ)||＋(１−α)||Ψ２Ｉ||１＋λ(ＸＴＤＸ)
ＸＴ：差分行列Ｘの転置
Ｄ：次の行列要素を有する対角行列のシステムノイズ行列
【００２９】
式（１１） Ｄｉｊ＝１／σｎ２ (ｉ＝ｊ)，０ (ｉ≠ｊ)
σ２：ノイズ分散及び本発明の実施で使用される画像装置でのノイズを表す変数
例えば、Ｘ線画像装置の場合、ノイズ変数σｎ２は、ｎｔｈのＸ線検出器に関連するノイズ分散である。
【００３０】
図１を参照して欲しい。本発明での方法の１つの実施で式（２）の目的関数を用いている。この場合、正則化変数αの初期化で開始している（ステップ１００）。正則化変数αの選択は、所望画像の散在と所望画像での先行画像の影響との間の取引を決定する。従って、正則化変数αの値は、目前の臨床用途に依存する。例えば、α≒０．３−０．７は、呼吸器系の画像用途には十分である。続いて、目的関数の式（２）の第１及び第２項が、各々初期化される（ステップ１０２、１０４）。第１項α||Ψ１(Ｉ−Ｉｐ)||の初期化は、ステップ１０６で始まり、先行画像Ｉｐは、所望画像Ｉの推定値から差し引いて、「差分画像」を生成する。先行画像Ｉｐの特定の選択及び所望画像Ｉの推定値は、画像化の手順及び特定の臨床用途に依存する。従って、これらの選択に対する異なる変更は、下記に詳細に説明する。差分画像は、次に、スパース変換Ψ１の適用で疎にさせる（ステップ１０８）。上記記載したように、スパース変換Ψ１は、ウェーブレット変換、１次差分、２次差分、及び、離散階調変換を含む何れの数学的演算となる。このスパース差分画像のＬ１ノルムは、ステップ１１０で計算される。この処理の結果は、ステップ１１２で表示される正則化変数αで重みを加える。
【００３１】
目的関数の式（２）の第２項(１−α)||Ψ２Ｉ||１の初期化は、ステップ１１４で開始し、所望画像Ｉの推測値は、スパース変換Ψ２の適用で疎にされる。続いて、このスパース所望画像推定値のＬ１ノルムがステップ１１６で計算される。離散階調変換∇ｍ，ｎがスパース変換Ψ２として選択されると、ステップ１１４、１１６で、所望画像推定値の全変数ＴＶの計算として見なされる。ＴＶは次式である。
式（１２） ＴＶ(Ｉ)＝||▽Ｉ||１＝Σ|▽Ｉ|
【００３２】
スパースした所望画像推定値のＬ１ノルムの計算後、結果を（１−α）で重みを付ける（ステップ１１８）。目的関数式（２）が、第１及び第２項を共に加えて形成される（ステップ１２０）。この目的関数は、例えば、非線形共役勾配法の使用で最小化される（ステップ１２２）。最小化処理は、停止基準が満たされるまで進められる。停止基準は、例えば、所望画像の現在推定値を前の繰返しからの所望画像の推定値との比較から成る。停止基準は、次式となる。
（式１３） ΣｉΣｊ(Ｉｉｊ（ｋ＋１）−Ｉｉｊ（ｋ）)２
Ｉｉｊ（ｋ＋１）：ピクセル位置（ｉ，ｊ）での所望画像の（ｋ＋１）ｔｈ推定値の値
Ｉｉｊ（ｋ）：ピクセル位置（ｉ，ｊ）での所望画像のｋｔｈ推定値の値
【００３３】
図２に示すように、本発明の方法の他の実施では目的関数式（７）を用い、正則化変数αの初期化で開始する（ステップ２００）。続いて、目的関数式（７）の第１及び第２項を各々初期化する（ステップ２０２、２０４）。この処理は、図１のステップ１０２、１０４と同様の方法で進行する。しかし、適合性条件式（１）は、ステップ２０６で初期化される第３項（λ||Ｘ||２２）に組み込まれている。最初に差分行列Ｘが生成される（ステップ２０８）。上記で詳細に記述したように、差分行列Ｘは、適合性条件式（１）に対応し、式（９）となる。
【００３４】
従って、差分行列Ｘは、システム行列Ａを所望画像Ｉの推定値に適用し、次に、所望画像に対応する取得された画像データＹを差し引いて確定される。差分行列ＸのＬ２ノルムの二乗が計算される（ステップ２１０）。差分行列ＸのＬ２ノルムの二乗が生成されると、ラグランジュ乗数λを決定し、差分行列Ｘに重みを付けるのに使用される（ステップ２１２）。上記のように、目前の臨床用途に応じて、ラグランジュ乗数の経験的な決定及び値が選択される。目的関数式（７）が、第１、第２及び第３項を共に加えて生成される（ステップ２２０）。目的関数は、非線形共役勾配法で最小化され（ステップ２２２）、最小化処理は、停止基準が満たされるまで進行する。
【００３５】
図３は、目的関数式（１０）を使用し、正則化変数αの初期化で開始する（ステップ３００）他の実施形態である。続いて、目的関数式（１０）の第１及び第２項が各々初期化される（ステップ３０２、３０４）。この処理は、図１に示すステップ１０２、１０４と同様の方法で進行する。しかし、適合性条件式（１）及び画像装置でのノイズの影響は、第３項λ(ＸＴＤＸ)に取り込んでいる（ステップ３０６）。最初に、差分行列Ｘが生成され（ステップ３０８）、図２のステップ２０８に記載するように、システムノイズ行列Ｄが生成される（ステップ３１０）。システムノイズ行列Ｄは、次式で決定される行列要素を有する対角行列である。
式（１４） Ｄｉｊ＝１／σｎ２ (ｉ＝ｊ)，０ (ｉ≠ｊ)
【００３６】
上記したように、σｎ２は、ノイズ分散で、本発明の実施で使用した画像装置のノイズを表す変数である。例えば、Ｘ線画像装置の場合、ノイズ変数σｎ２は、ｎｔｈのＸ線検出器に関連するノイズ分散である。システムノイズ行列Ｄが生成されると、次の行列乗算が実行される（ステップ３１２）。
（式１５） ＸＴＤＸ
【００３７】
この演算の結果は、続いて、ラグランジュ乗数で重みが付けられる（ステップ３１４）。目的関数式（１０）は、続いて、第１、第２及び第３項を共に加えて形成される（ステップ３２０）。目的関数は、例えば非線形共役勾配法の使用で最小化される（ステップ３２２）。この最小化処理は、停止基準が満たされるまで進行する。
【００３８】
ＣアームＸ線画像装置
【００３９】
本発明の実施形態の画像再構成に使用したＸ線装置（特に、図４Ａ及び図４Ｂ）は、治療処置に関連した使用に設計されている。一端部でＸ線源装置４１２と他端部でＸ線検出器配列装置４１４とを担持するＣアーム４１０を有する架台が特徴である。架台は、台４１６上の患者の回りで、Ｘ線源４１２及び検出器４１４を異なる位置及び角度に向けることが可能で、医師の患者への接近を可能にしている。
【００４０】
架台は、台４１６下に延びている水平脚４２０と、台４１６から離れ、水平脚４２０の端部で上部に延在する垂直脚４２２と、を有しているＬ形状の台座４１８を有している。 支持腕４２４は、垂直脚４２２の上端に回転可能に取付けられ、水平の旋回軸４２６に対し回転できる。旋回軸４２６は、台４１６の中心線に並んでおり、腕４２４は、旋回軸４２６から径方向上部に延在し、その外端部でＣアーム駆動装置４２７を支持している。Ｃアーム４１０は、駆動装置４２７に摺動可能に固定され、駆動モータ（図示せず）に結合されている。駆動モータは、Ｃアーム４１０を摺動させ、Ｃ軸４２８に対し回転させる（矢印４３０）。旋回軸４２６及びＣ軸４２８は、台４１６上に位置する等角点で互いに垂直に交差している。
【００４１】
Ｘ線源装置４１２は、Ｃアーム４１０の一端部に取付けられ、検出器配列装置４１４は、他端部に取付けられている。下記に詳細に記載するが、Ｘ線源４１２は、検出器４１４に向かう円錐状ビームのＸ線を放射する。両装置４１２、４１４は、旋回軸４２６に径方向内側に延在し、円錐状ビームの中心光が装置の等角点４３６を通過する。円錐状ビームの中心光を、装置の等角点に対し、台４１６上に載せた被検者からのＸ線減衰データ取得中、旋回軸４２６又はＣ軸４２８の周り、又は、双方で回転できる。
【００４２】
図５Ａに示すように、Ｘ線源装置４１２は、電圧を加えると円錐状４３３のＸ線ビームを放射するＸ線源４３２を備えている。中心光４３４は、装置の等角点４３６を通過し、検出装置４１４に格納されている２次元平坦パネルのデジタル検出器４３８に当たる。検出器４３８は、２０４８×２０４８要素の２次元配列の検出器素子で、４１ｃｍ×４１ｃｍの寸法を有している。各素子は、患者を通過して突き当たったＸ線強度及びＸ線の減衰を表す電気信号を生成する。走査中、Ｘ線源４３２及び検出器配列４３８は、装置の等角点４３６で回転し、異なる角度からの減衰投影データを取得する。検出器配列は、毎秒３０の投影又は視野の取得が可能で、これが幾つの視野が所定の走査通過と速度で取得される制限要因となる。
【００４３】
装置４１２、４１４の回転及びＸ線源４３２の操作は、ＣＴ装置の制御機構部４４０で制御されている（図４Ｂ参照）。制御機構部４４０は、Ｘ線源４３２に電力及びタイミング信号を与えるＸ線制御部４４２を含んでいる。制御機構部４４０のデータ取得装置（ＤＡＳ）４４４は、検出器素子４３８からのデータを採取し、そのデータを画像再構成部４４５に送る。画像構成部４４５は、ＤＳＡ４４４からデジタル化したＸ線データを受け取り、本発明の方法で高速の画像再構成を実行する。再構成された画像は、大容量記憶デバイス４４９に画像を蓄積又は画像を更に処理するコンピュータ４４６に入力として適用される。
【００４４】
制御機構部４４０は、旋回軸モータ制御部４４７及びＣ軸モータ制御部４４８も含んでいる。コンピュータ４４６からの動作コマンドに応じて、モータ制御部４４７、４４８は、旋回軸４２６及びＣ軸４２８で各々の回転を生成するＸ線装置のモータに電力を供給する。コンピュータ４４６で実行されるプログラムは、モータ駆動部４４７に動作コマンドを生成し、所定の走査経路で装置４１２、４１４を動かす。
【００４５】
コンピュータ４４６は、キーボード及び他の手動操作制御部を有するコンソール４５０を介した操作者からのコマンド及び走査変数を受け取る。関与する陰極線管表示部は、操作者の再構成された画像及びコンピュータ４４６からのほかのデータの観察が可能である。操作者が与えたコマンドは、蓄積されたプログラムの指示に従うコンピュータ４４６で使用され、ＤＡＳ４４４、Ｘ線制御部４４２及びモータ制御部４４８に制御信号及び情報を与える。加えて、コンピュータ４４６は、装置の等角点４３６に患者を位置付ける電動台４１６を制御する台のモータ制御部４５４を操作する。
【００４６】
上記記載のＸ線装置は、蛍光透視法で操作し実時間で２次元画像を生成できる。主として、このモードの使用は、患者にカテーテル法を実施し、蛍光透視法を処置の誘導に使用する場合である。誘導されている器具及びその周囲組織が十分に見えるＳＮＲの画像フレームが生成されるのが好ましい。更に、実時間に近い処置画像を与えるフレーム速度及び可能な限り低いＸ線量となる画像フレームの形成が望ましい。
【００４７】
磁気共鳴画像装置
【００４８】
図６は、本発明の実施形態のＭＲＩ装置を示す。ＭＲＩ装置は、表示部６１２及びキーボード６１２を有するワークステーション６１０を含む。ワークステーション６１０は、市販で入手可能なＯＳで動作する市販で入手可能な機械の演算部６１６を含む。ワークステーション６１０は、ＭＲＩ装置に走査指示の登録を可能とする操作インターフェイスを与える。ワークステーション６１０は、パルス系列サーバ６１８、データ取得サーバ６２０、データ処理サーバ６２２、及び、データ蓄積サーバ６２３の４つのサーバに接続されている。ワークステーション６１０及び各サーバ６１８、６２０、６２２、６２３は、互いに通信可能に接続されている。
【００４９】
パルス系列サーバ６１８は、ワークステーション６１０からダウンロードした命令に応じて作用し、傾斜装置６２４及びＲＦ装置６２６を操作する。所定の走査の実行に必要な傾斜波形が形成され、傾斜装置６２４に適用される。傾斜装置６２４は、装置６２８内の傾斜コイルを励起し、磁場傾斜Ｇｘ、Ｇｙ、Ｇｚを生成し、ＭＲ信号の位置エンコードに使用される。傾斜コイル装置６２８は、分極磁石６３２及びＲＦコイル６３４全体を含む磁石装置６３０の部分を形成する。
【００５０】
ＲＦ励起波形は、ＲＦ装置６２６によりＲＦコイル６３４に加えられ、所定の磁気共鳴パルス系列を実施する。ＲＦコイル６３４又は別個の局在コイル（図６に図示せず）で検出される応答ＭＲ信号は、パルス系列サーバ６１８で形成されるコマンドの指令で、ＲＦ装置６２６での受け取り、増幅、復調、フィルター及びデジタル化される。ＲＦ装置６２６は、ＭＲパルス系列に使用される幅広いＲＦパルスを形成するＲＦ送信器を含む。ＲＦ送信器は、パルス系列サーバ６１８からの走査指示及び走査指令に応じ、所望の周波数、位相及びパルス振幅波形を形成する。生成されたＲＦパルスは、ＲＦコイル全体、又は、１つ又はそれ以上の局在コイル又はコイル配列（図６に図示せず）に加えられる。
【００５１】
ＲＦ装置６２６は、１つ又はそれ以上のＲＦ受信路も含んでいる。各ＲＦ受信路は、接続されたコイルで受信されるＭＲ信号を増幅するＲＦ増幅器と、受信ＭＲ信号のＩ及びＱ直交成分を検出及びデジタル化する検出器とを含んでいる。受信されたＭＲ信号の大きさは、採取点でのＩ及びＱ成分の二乗の和の平方根、
Ｍ＝√(Ｉ２＋Ｑ２)
及び、受信ＭＲ信号の位相で決定される。
φ＝ｔａｎ−１(Ｑ/Ｉ)
【００５２】
パルス系列サーバ６１８は、生理学的取得制御部６３６から患者データを適宜受け取る。制御部６３６は、患者に接続された幾つかの異なるセンサ（例、電極からのＥＣＧ信号又はベローズからの呼吸信号）からの信号を受け取る。これらの信号は、パルス系列サーバ６１８で使用され、被検者の呼吸又は心臓の鼓動で走査の実行を同期又はゲートさせる。
【００５３】
パルス系列サーバ６１８は、走査室インターフェイス回路６３８にも接続している。インターフェイス回路６３８は、患者の状態及び磁気装置に関連する種々のセンサから信号を受け取る。走査室インターフェイス回路６３８を経て、患者位置決め装置６４０は、走査中、患者を所望の位置に動かすコマンドを受け取る。
【００５４】
ＲＦ装置６２６で形成されたデジタルＭＲ信号サンプルは、データ取得サーバ６２０で受信される。データ取得サーバ６２０は、ワークステーション６１０からダウンロードされる指示に呼応して、実時間ＭＲデータを受け取り、過剰データによるデータ損失を防ぐバッファ蓄積部を提供する。幾つかの走査で、データ取得サーバ６２０は、取得したＭＲデータをデータ処理サーバ６２２に送るだけの場合がある。しかし、取得したＭＲデータから誘導される情報で、走査の更なる実行を制御する走査の場合、データ取得サーバ６２０は、そのような情報を形成するようにプログラムされ、情報をパルス系列サーバ６１８に運ぶ。例えば、予備走査中、ＭＲデータが取得され、そのＭＲデータは、パルス系列サーバ６１８で実行されるパルス系列の較正に使用される。又、自動信号を走査中に取得させ、ＲＦ又は傾斜装置の操作変数の調整、又は、ｋ空間が採取される視野順の制御に使用できる。データ取得サーバ６２０は、磁気共鳴血管造影（ＭＲＡ）走査での造影剤の到達の検出に使用されるＭＲ信号の処理に使用しても良い。全ての実施例で、データ取得サーバ６２０は、ＭＲデータを取得し、そのデータを実時間で処理し、走査の制御に使用される情報を形成する。
【００５５】
データ処理サーバ６２２は、データ取得サーバ６２０からＭＲデータを受け取り、ワークステーション６１０からダウンロードされた指示に従い、そのデータを処理する。処理は、例えば、２又は３次元画像を形成するための未処理ｋ空間ＭＲデータのフーリエ変換、フィルタの再構成画像への適用、取得ＭＲデータの逆投影法画像再構成の実施、機能ＭＲ画像の計算、及び、動画又は流れる画像の計算等を含む。
【００５６】
データ処理サーバ６２２で再構成された画像は、ワークステーション６１０に戻されて蓄積される。実時間画像が、データベースメモリキャッシュ（図示せず）に蓄積される。キャッシュから、磁石装置６３０近傍に位置する操作表示部６１２又は表示部６４２に向けて、担当医の使用に出力される。バッチモード画像又は選択された実時間画像は、ディスク蓄積部６４４上のホストデータベースに蓄積される。画像が再構成され、蓄積部に転送されると、データ処理サーバ６２２は、ワークステーション６１０のデータ蓄積サーバ６２３に通知する。操作者は、画像の保管、フィルムの形成、又は、画像をネットワーク経由で他の施設に送るためにワークステーション６１０を使用する。
【００５７】
従来のＭＲＩ装置を本発明で使用できるが、１つの実施形態では、医者の接近を可能にする設計にしたＭＲＩ装置を使用する。特に図７に示すように、術中ＭＲ画像処理が実施される時、患者は、磁石装置６３０内に配置され、患者の関与領域を、２つの離れた磁石の輪６９０、６９２間に位置する装置の等角点近傍に位置させる。磁石の輪６９０、６９２間に立つ医師は、患者の重要な領域への接近に制限されない。処置中、患者は動いても良く、医師が動かしても良い。集束超音波（ＦＵＳ）切除装置６８０に接続された超音波振動子（図示せず）を医師が操作し、患者の選択した領域を治療しても良い。表示部６４２に形成される画像を、医師が、切除装置の支援及び治療終了の決定に使用できる。
【００５８】
画像再構成
【００５９】
図８を参照されたい。本発明の実施形態は、画像装置で画像データの取得で開始する（ステップ８００）。図４Ａ、４Ｂに示すようなＣアームＸ線画像装置で実施する場合、画像データは、Ｘ線投影空間を採取する一連のＸ線投影として取得される。図６、７に示すようなＭＲＩ装置で実施する場合、画像データは、ｋ空間を採取する一連のｋ空間投影として取得される。１つの実施形態として、取得された画像データは、時間分解画像データで、画像フレームの時系列に対応する。より具体的には、「画像フレーム」は、対応する時点で取得される画像データ又は「タイムフレーム」から再構成される。
【００６０】
取得された画像データから、次に、初期画像フレームが、各々の対応するタイムフレームに対し再構成される（ステップ８０２）。これは標準的な画像再構成である。例えば、画像データがＸ線画像データの場合、フィルタ補正逆投影の画像再構成法が使用される。
例えば、画像データが一連のｋ空間投影の場合、逆投影法も使用される。他の方法として、ｋ空間投影データを直交座標格子に再編し、従来のフーリエ変換再構成法の使用が出来る。各タイムフレームで取得される画像データは、下取りの採取で、画像の縞、画像のギザギザ等の画像の不自然さが、上記の再構成された初期画像フレームに存在する。これら初期画像フレームは、被検者の先行画像の形成に使用される（ステップ８０４）。例えば、先行画像は、入手できる初期画像フレームの平均化で形成される。
【００６１】
続いて、得られた先行画像、取得画像データ、及び、現在入手できる初期画像フレームを使用して、対応する画像フレームの推定値を形成する（ステップ８０６）。画像フレーム推定値を、例えば、幾何学的再構成技術（ＡＲＴ）を使用して形成する。これらの画像フレーム推定値を利用して、先行画像と共に、上記記載の最小化法で被検者の最終画像フレームを再構成する（ステップ８０８）。特に、次の目的関数が反復して最小化される。
||Ψ１(Ｉａｖｅ−Iｐ)||１＋α||Ψ２Ｉｔ||１＋βΣｔ＝２Ｎ||Ψ３(Iｔ−Iｔ−１)||１
Ψ１、Ψ２、Ψ３：スパース変換
Ｉａｖｅ：現在入手可能画像フレーム推定値Ｉｔ
Ｉｐ：先行画像
α、β：正則化変数
Ｉｔ：ｔｔｈの画像フレーム
上覧の目的関数は、式（２）、式（７）、式（１０）と少なくとも２点で異なる。第１に、目的関数の最後の項は、
βΣｔ＝２Ｎ||Ψ３(Ｉｔ−Ｉｔ−１)||１
入手可能な画像フレーム推定値の全てに及ぶ２つの隣接する画像フレーム推定値間の差異の和である。このように、現在の画像フレームＩｔに対する最近接画像フレームＩｔ−１は、
式（２）、式（７）、式（１０）の先行画像と同様に作用する。下取りデータから現在の画像フレームの再構成を制約するために使用される。第２として、目的関数の第１項は、
||Ψ１(Ｉａｖｅ−Iｐ)||１
先行画像と入手可能な画像フレーム推定値の平均との差異を計算する。この理由から、前述の目的関数を用いる本発明の実施形態は、実時間画像の適用に適している。実時間画像の場合、データは連続して取得され、データの取得中に画像が再構成される。この方法は、例えば、Ｘ線又はＭＲ画像装置で誘導される治療処置中で取得される画像の品質の向上を可能にする。
【００６２】
適切な停止基準に達すると（ステップ８１０）、画像再構成法は中止され、所望の画像フレームが蓄積又はモニタに表示される。停止基準が満たされない場合、画像フレーム推定値は、目的関数の最小化の出力で更新される（ステップ８１２）。続いて、これら更新された画像フレーム推定値は、上覧の目的関数で利用される。典型的な停止基準は、２つの隣接する繰返しの差異が前以て設定された閾値よりも小さい場合である。
【００６３】
本発明は、１つ又はそれ以上の実施形態で記載している。これら記載したものの等価物、変更、改変は可能で、本発明の範囲内になる。
【技術分野】
【０００１】
本願は、表題"Method for Image Reconstruction Using Prior Image Constrained Compressed Sensing"で出願(2008年1月14日付)している米国仮出願(番号61/020,847)、及び、表題"Method for Image Reconstruction Using Prior Image Constrained Compressed Sensing"で(2008年6月9日付)出願している米国仮出願(番号61/059,891)において優先権を主張している。
【０００２】
この発明は、米国政府機関の国立衛生研究所の支援(NIH EB005712及びNIH EB007021)でなされた。米国政府は、本発明に特定の権利を有する。
【背景技術】
【０００３】
発明の領域は、医用画像に関し、特に画像を再構成する方法に関する。
【０００４】
コンピュータ断層撮影装置において、Ｘ線源は、直交座標系のＸ−Ｙ面（画像面）内で平行になる円錐形状ビームを放射する。Ｘ線ビームは、撮像する被写体（例、内科患者）を通過し、配列した放射線検出器に突き当たる。透過放射線強度は、被写体のＸ線ビームの減衰に依存し、各検出器は、ビーム減衰の測定となる個別の電気信号を生ずる。全ての検出器からの減衰測定は、個別に取得され、透過像、減衰像、又は、投影と呼ばれるものを形成する。
【０００５】
従来のＣＴ装置の光源及び検出器配列は、架台上で画像面内及び被写体周囲で回転させ、 Ｘ線ビームの被写体に交差する角度を常に変化させる。任意の角度での検出器配列からの透過像は、１つの視野となり、被写体の走査は、Ｘ線源及び検出器の一回転での異なる角度方向で形成される一組の視野を含む。２Ｄ走査の場合、データは、被写体から得られる２次元薄片に対応する画像の構成に処理される。２Ｄデータから画像を再構成する一般的な方法は、フィルタ補正逆投影法と称されている。この画像再構成処理は、走査中に取得される減衰測定を、表示部の対応するピクセルの輝度制御に使用されるＣＴ数又はＣＴ値（Hounsfield units）と呼ばれる整数値に変換される。
【発明の概要】
【発明が解決しようとする課題】
【０００６】
本発明は、画像再構成の方法を提供するものである。
【課題を解決するための手段】
【０００７】
本発明の画像再構成法は、正確な再構成の利点を信号対ノイズ比（ＳＮＲ）のゲインと組合せている。つまり、被検者の適切な先行画像を画像再構成処理に組み込んでいる。ＳＮＲの増加に加えて、本発明の方法は、分解能を一次的に増す。更に、得られる所望画像及び投影の数に対し、本発明の方法は、従来の画像再構成法（例、圧縮センシング：ＣＳ）が失敗している正確な画像再構成を与える。
【０００８】
本発明は、時間分解画像データから良質の画像フレームを実時間で再構成する方法を提供する。この方法は、実時間画像方法を利用する画像誘導治療処置に適用できる。例えば、Ｘ線蛍光誘導治療処置は、本発明の方法で利点を得る。この場合、良質の画像フレームが、患者への全Ｘ線量の減少で再構成ができる。
【０００９】
本発明の方法は、時系列の画像フレームでの特定画像フレームの再構成を与える。本方法は、時系列の最も近い画像フレームから現在の画像フレームを差し引く。例えば、前画像フレームが、現在の画像フレームから差し引かれる。このように、各新規に再構成された画像フレームが、動的な先行画像として使用され、現在の画像フレームを疎にする。現在の画像フレームは、更に少ない投影データから正確な再構成を可能にする。
【００１０】
本発明の利点を下記明細書で説明する。記載において、添付図を用いて発明の実施形態を説明している。実施形態は、発明の範囲を限定するものでなく、発明の解釈の参照とするものである。
【図面の簡単な説明】
【００１１】
【図１】本発明で利用している画像再構成法の実施形態の工程図である。
【図２】本発明で利用している画像再構成法の他の実施形態の工程図である。
【図３】本発明で利用している画像再構成法の別の実施形態の工程図である。
【図４】（Ａ）本発明を使用するＣアームＸ線装置の図。（Ｂ）図４ＡのＣアームＸ線装置の構成図である。
【図５】（Ａ）図４ＡのＣアームＸ線装置のＸ線源及び検出器の図。（Ｂ）図４ＡのＣアームＸ線装置で用いられるＣアーム走査経路の図である。
【図６】図６は、本発明を使用する磁気共鳴画像（ＭＲＩ）装置の構成図である。
【図７】図７は、本発明を使用する他の磁気共鳴画像（ＭＲＩ）装置の構成図である。
【図８】図８は、本発明の実施形態の工程を説明する工程図である。
【発明を実施するための形態】
【００１２】
一般的に、一組のデータから画像を再構成する方法は、測定データサンプルＹから所望の画像Ｉを推測する一連の計算ステップを含む。更に具体的には、画像再構成は、次の適合性条件を満たす。
【００１３】
式（１） ＡＩ＝Ｙ
【００１４】
Ａはシステム行列で、一般的には、所望画像Ｉを取得データサンプルＹに関連付ける前方投影演算子として見なせる。コンピュータ断層撮影（ＣＴ）画像の処理の場合、システム行列は再投影操作を含む。式（１）の適合性条件は、別の言葉で言えば、画像が忠実に再構成される時、前方演算は、測定した投影データの正確な推定値を生成するために、実データ取得手続に略似ているべきである。
【００１５】
本発明の実施で使用する画像再構成法は、良質の画像を提供する。「先行画像」は、圧縮センシング（ＣＳ）原理が利用されている反復画像再構成法を抑制するために使用されている。例えば、ＣＳで一般に使用されているスパース変換に加えて、画像を所望画像から先行画像を取り去って更に疎にする。その結果、ＣＳ法で要求されるよりも少ないサンプル数の使用で、画像が正確に再構成される。
【００１６】
更に具体的には、先行画像Ｉｐ及び再構成する所望画像Ｉにおいて、画像再構成の本発明の方法は、次の目的関数の最小化で達成される。
【００１７】
式（２） α||Ψ１(Ｉ−Ｉｐ)||１＋(１−α)||Ψ２Ｉ||1
Ψ１、Ψ２：スパース変換
||...||１：Ｌ１ノルム演算
α：目的関数式（２）の２つの項の相対的重みの制御に利用される正則化変数
【００１８】
上記から、次式
式（３） ||x||１＝Σｉ＝１Ｎ|xｉ|
は、Ｎ次元ベクトルｘのＬ１ノルムを意味する。一般的には、真のＬ１ノルムからの偏差は、所望画像に適切な画質を維持する限り可能である。例えば、目的関数式（２）は、次式に一般化できる。
【００１９】
式（４） α||Ψ１(Ｉ−Ｉｐ)||ｐ＋(１−α)||Ψ２Ｉ||ｐ
||...||ｐ：次の形式を有するＬｐノルム演算
式（５） ||ｘ||ｐ＝(Σｉ＝１Ｎ|ｘｊ|ｐ)１/ｐ
【００２０】
上記のように、好ましくはｐ＝１．０であるが、しかし、異なるｐ値が可能である。
【００２１】
当業者には周知であるが、ｐ値が１．０から更に外れると、再構成した所望画像の劣化が一層明白になる。
【００２２】
式（２）でのスパース変換Ψ１、Ψ２は一般的には異なるが、同一であっても良い。例次的なスパース変換は、ウェーブレット変換、１次差分、２次差分、及び、次式の離散階調変換を含んでいる。
式（６） ▽m,nＩ(m,n)＝√([I(m+1,n)−I(m,n)]２＋[I(m,n＋1)−Ｉ(m,n)]２)
指数ｍ、ｎ：画像Ｉのピクセルの位置
▽m,nＩ(m,n)で特定される画像は、通常、「階調画像」と称されている。
【００２３】
目的関数式（２）の両方の項が重要である。重要性の結果として、正則化変数の選択αが、全体の画像再構成処理の制御に利用される。正則化変数αの選択は、先行画像Ｉｐの選択及び手元の医療用途の選択にも依存する。例えば、目的関数式（２）の第２項は、先行画像Ｉｐから潜在的に受け継ぐ縞のノイズを軽減する。実例として、α≒０．３−０．７の正則化変数の選択は、呼吸器官の画像適用には一般的に十分である。
【００２４】
適合性条件式（１）を全体の画像再構成に一層良く取り込むために、ラグランジュ乗数法が利用される。この場合、適合性条件は、式（２）で設定した目的関数の最小化に更に制約を加えて使用される。新しい目的関数は、次式の形式となる。
式（７） α||Ψ１(Ｉ−Ｉｐ)||1＋(１−α)||Ψ２Ｉ||1＋λ||Ｘ||２２
λ：ラグランジュ乗数
Ｘ：差分行列
||...||２２：平方Ｌ２ノルム演算
Ｎ次元ベクトルｘに対し、Ｌ２ノルムは次式となる。
式（８） ||ｘ||２２＝Σｉ＝１Ｎｘｉ２
【００２５】
式（７）の差分行列は、適合性条件式（１）から成り、次式となる。
式（９） Ｘ＝ＡＩ−Ｙ
【００２６】
ラグランジュ乗数λは、本発明の実施で使用される特定の画像装置に対し経験的に決定される。例えば、ラグランジュ乗数λは、所望のデータ適合性条件と先行画像Ｉｐの類似性との間で所定の交換条件で決定される。大きいラグランジュ乗数λが選択されると、再構成された画像は、低ノイズ分散を有する。しかし、先行画像に特徴な高い空間分解能の損失となる。小さいラグランジュ乗数λが使用されると、先行画像の特徴である高い空間分解能が維持されるが、ノイズ分散が所望の画像で高くなる。このような状況は、利用される画像装置で達成されるコントラスト対ノイズ比に影響を及ぼす。
【００２７】
式（７）に示す目的関数は、画像装置のノイズを説明するために変更できる。この場合、次の目的関数が最小化される。
【００２８】
式（１０） α||Ψ１(Ｉ−Ｉｐ)||＋(１−α)||Ψ２Ｉ||１＋λ(ＸＴＤＸ)
ＸＴ：差分行列Ｘの転置
Ｄ：次の行列要素を有する対角行列のシステムノイズ行列
【００２９】
式（１１） Ｄｉｊ＝１／σｎ２ (ｉ＝ｊ)，０ (ｉ≠ｊ)
σ２：ノイズ分散及び本発明の実施で使用される画像装置でのノイズを表す変数
例えば、Ｘ線画像装置の場合、ノイズ変数σｎ２は、ｎｔｈのＸ線検出器に関連するノイズ分散である。
【００３０】
図１を参照して欲しい。本発明での方法の１つの実施で式（２）の目的関数を用いている。この場合、正則化変数αの初期化で開始している（ステップ１００）。正則化変数αの選択は、所望画像の散在と所望画像での先行画像の影響との間の取引を決定する。従って、正則化変数αの値は、目前の臨床用途に依存する。例えば、α≒０．３−０．７は、呼吸器系の画像用途には十分である。続いて、目的関数の式（２）の第１及び第２項が、各々初期化される（ステップ１０２、１０４）。第１項α||Ψ１(Ｉ−Ｉｐ)||の初期化は、ステップ１０６で始まり、先行画像Ｉｐは、所望画像Ｉの推定値から差し引いて、「差分画像」を生成する。先行画像Ｉｐの特定の選択及び所望画像Ｉの推定値は、画像化の手順及び特定の臨床用途に依存する。従って、これらの選択に対する異なる変更は、下記に詳細に説明する。差分画像は、次に、スパース変換Ψ１の適用で疎にさせる（ステップ１０８）。上記記載したように、スパース変換Ψ１は、ウェーブレット変換、１次差分、２次差分、及び、離散階調変換を含む何れの数学的演算となる。このスパース差分画像のＬ１ノルムは、ステップ１１０で計算される。この処理の結果は、ステップ１１２で表示される正則化変数αで重みを加える。
【００３１】
目的関数の式（２）の第２項(１−α)||Ψ２Ｉ||１の初期化は、ステップ１１４で開始し、所望画像Ｉの推測値は、スパース変換Ψ２の適用で疎にされる。続いて、このスパース所望画像推定値のＬ１ノルムがステップ１１６で計算される。離散階調変換∇ｍ，ｎがスパース変換Ψ２として選択されると、ステップ１１４、１１６で、所望画像推定値の全変数ＴＶの計算として見なされる。ＴＶは次式である。
式（１２） ＴＶ(Ｉ)＝||▽Ｉ||１＝Σ|▽Ｉ|
【００３２】
スパースした所望画像推定値のＬ１ノルムの計算後、結果を（１−α）で重みを付ける（ステップ１１８）。目的関数式（２）が、第１及び第２項を共に加えて形成される（ステップ１２０）。この目的関数は、例えば、非線形共役勾配法の使用で最小化される（ステップ１２２）。最小化処理は、停止基準が満たされるまで進められる。停止基準は、例えば、所望画像の現在推定値を前の繰返しからの所望画像の推定値との比較から成る。停止基準は、次式となる。
（式１３） ΣｉΣｊ(Ｉｉｊ（ｋ＋１）−Ｉｉｊ（ｋ）)２
Ｉｉｊ（ｋ＋１）：ピクセル位置（ｉ，ｊ）での所望画像の（ｋ＋１）ｔｈ推定値の値
Ｉｉｊ（ｋ）：ピクセル位置（ｉ，ｊ）での所望画像のｋｔｈ推定値の値
【００３３】
図２に示すように、本発明の方法の他の実施では目的関数式（７）を用い、正則化変数αの初期化で開始する（ステップ２００）。続いて、目的関数式（７）の第１及び第２項を各々初期化する（ステップ２０２、２０４）。この処理は、図１のステップ１０２、１０４と同様の方法で進行する。しかし、適合性条件式（１）は、ステップ２０６で初期化される第３項（λ||Ｘ||２２）に組み込まれている。最初に差分行列Ｘが生成される（ステップ２０８）。上記で詳細に記述したように、差分行列Ｘは、適合性条件式（１）に対応し、式（９）となる。
【００３４】
従って、差分行列Ｘは、システム行列Ａを所望画像Ｉの推定値に適用し、次に、所望画像に対応する取得された画像データＹを差し引いて確定される。差分行列ＸのＬ２ノルムの二乗が計算される（ステップ２１０）。差分行列ＸのＬ２ノルムの二乗が生成されると、ラグランジュ乗数λを決定し、差分行列Ｘに重みを付けるのに使用される（ステップ２１２）。上記のように、目前の臨床用途に応じて、ラグランジュ乗数の経験的な決定及び値が選択される。目的関数式（７）が、第１、第２及び第３項を共に加えて生成される（ステップ２２０）。目的関数は、非線形共役勾配法で最小化され（ステップ２２２）、最小化処理は、停止基準が満たされるまで進行する。
【００３５】
図３は、目的関数式（１０）を使用し、正則化変数αの初期化で開始する（ステップ３００）他の実施形態である。続いて、目的関数式（１０）の第１及び第２項が各々初期化される（ステップ３０２、３０４）。この処理は、図１に示すステップ１０２、１０４と同様の方法で進行する。しかし、適合性条件式（１）及び画像装置でのノイズの影響は、第３項λ(ＸＴＤＸ)に取り込んでいる（ステップ３０６）。最初に、差分行列Ｘが生成され（ステップ３０８）、図２のステップ２０８に記載するように、システムノイズ行列Ｄが生成される（ステップ３１０）。システムノイズ行列Ｄは、次式で決定される行列要素を有する対角行列である。
式（１４） Ｄｉｊ＝１／σｎ２ (ｉ＝ｊ)，０ (ｉ≠ｊ)
【００３６】
上記したように、σｎ２は、ノイズ分散で、本発明の実施で使用した画像装置のノイズを表す変数である。例えば、Ｘ線画像装置の場合、ノイズ変数σｎ２は、ｎｔｈのＸ線検出器に関連するノイズ分散である。システムノイズ行列Ｄが生成されると、次の行列乗算が実行される（ステップ３１２）。
（式１５） ＸＴＤＸ
【００３７】
この演算の結果は、続いて、ラグランジュ乗数で重みが付けられる（ステップ３１４）。目的関数式（１０）は、続いて、第１、第２及び第３項を共に加えて形成される（ステップ３２０）。目的関数は、例えば非線形共役勾配法の使用で最小化される（ステップ３２２）。この最小化処理は、停止基準が満たされるまで進行する。
【００３８】
ＣアームＸ線画像装置
【００３９】
本発明の実施形態の画像再構成に使用したＸ線装置（特に、図４Ａ及び図４Ｂ）は、治療処置に関連した使用に設計されている。一端部でＸ線源装置４１２と他端部でＸ線検出器配列装置４１４とを担持するＣアーム４１０を有する架台が特徴である。架台は、台４１６上の患者の回りで、Ｘ線源４１２及び検出器４１４を異なる位置及び角度に向けることが可能で、医師の患者への接近を可能にしている。
【００４０】
架台は、台４１６下に延びている水平脚４２０と、台４１６から離れ、水平脚４２０の端部で上部に延在する垂直脚４２２と、を有しているＬ形状の台座４１８を有している。 支持腕４２４は、垂直脚４２２の上端に回転可能に取付けられ、水平の旋回軸４２６に対し回転できる。旋回軸４２６は、台４１６の中心線に並んでおり、腕４２４は、旋回軸４２６から径方向上部に延在し、その外端部でＣアーム駆動装置４２７を支持している。Ｃアーム４１０は、駆動装置４２７に摺動可能に固定され、駆動モータ（図示せず）に結合されている。駆動モータは、Ｃアーム４１０を摺動させ、Ｃ軸４２８に対し回転させる（矢印４３０）。旋回軸４２６及びＣ軸４２８は、台４１６上に位置する等角点で互いに垂直に交差している。
【００４１】
Ｘ線源装置４１２は、Ｃアーム４１０の一端部に取付けられ、検出器配列装置４１４は、他端部に取付けられている。下記に詳細に記載するが、Ｘ線源４１２は、検出器４１４に向かう円錐状ビームのＸ線を放射する。両装置４１２、４１４は、旋回軸４２６に径方向内側に延在し、円錐状ビームの中心光が装置の等角点４３６を通過する。円錐状ビームの中心光を、装置の等角点に対し、台４１６上に載せた被検者からのＸ線減衰データ取得中、旋回軸４２６又はＣ軸４２８の周り、又は、双方で回転できる。
【００４２】
図５Ａに示すように、Ｘ線源装置４１２は、電圧を加えると円錐状４３３のＸ線ビームを放射するＸ線源４３２を備えている。中心光４３４は、装置の等角点４３６を通過し、検出装置４１４に格納されている２次元平坦パネルのデジタル検出器４３８に当たる。検出器４３８は、２０４８×２０４８要素の２次元配列の検出器素子で、４１ｃｍ×４１ｃｍの寸法を有している。各素子は、患者を通過して突き当たったＸ線強度及びＸ線の減衰を表す電気信号を生成する。走査中、Ｘ線源４３２及び検出器配列４３８は、装置の等角点４３６で回転し、異なる角度からの減衰投影データを取得する。検出器配列は、毎秒３０の投影又は視野の取得が可能で、これが幾つの視野が所定の走査通過と速度で取得される制限要因となる。
【００４３】
装置４１２、４１４の回転及びＸ線源４３２の操作は、ＣＴ装置の制御機構部４４０で制御されている（図４Ｂ参照）。制御機構部４４０は、Ｘ線源４３２に電力及びタイミング信号を与えるＸ線制御部４４２を含んでいる。制御機構部４４０のデータ取得装置（ＤＡＳ）４４４は、検出器素子４３８からのデータを採取し、そのデータを画像再構成部４４５に送る。画像構成部４４５は、ＤＳＡ４４４からデジタル化したＸ線データを受け取り、本発明の方法で高速の画像再構成を実行する。再構成された画像は、大容量記憶デバイス４４９に画像を蓄積又は画像を更に処理するコンピュータ４４６に入力として適用される。
【００４４】
制御機構部４４０は、旋回軸モータ制御部４４７及びＣ軸モータ制御部４４８も含んでいる。コンピュータ４４６からの動作コマンドに応じて、モータ制御部４４７、４４８は、旋回軸４２６及びＣ軸４２８で各々の回転を生成するＸ線装置のモータに電力を供給する。コンピュータ４４６で実行されるプログラムは、モータ駆動部４４７に動作コマンドを生成し、所定の走査経路で装置４１２、４１４を動かす。
【００４５】
コンピュータ４４６は、キーボード及び他の手動操作制御部を有するコンソール４５０を介した操作者からのコマンド及び走査変数を受け取る。関与する陰極線管表示部は、操作者の再構成された画像及びコンピュータ４４６からのほかのデータの観察が可能である。操作者が与えたコマンドは、蓄積されたプログラムの指示に従うコンピュータ４４６で使用され、ＤＡＳ４４４、Ｘ線制御部４４２及びモータ制御部４４８に制御信号及び情報を与える。加えて、コンピュータ４４６は、装置の等角点４３６に患者を位置付ける電動台４１６を制御する台のモータ制御部４５４を操作する。
【００４６】
上記記載のＸ線装置は、蛍光透視法で操作し実時間で２次元画像を生成できる。主として、このモードの使用は、患者にカテーテル法を実施し、蛍光透視法を処置の誘導に使用する場合である。誘導されている器具及びその周囲組織が十分に見えるＳＮＲの画像フレームが生成されるのが好ましい。更に、実時間に近い処置画像を与えるフレーム速度及び可能な限り低いＸ線量となる画像フレームの形成が望ましい。
【００４７】
磁気共鳴画像装置
【００４８】
図６は、本発明の実施形態のＭＲＩ装置を示す。ＭＲＩ装置は、表示部６１２及びキーボード６１２を有するワークステーション６１０を含む。ワークステーション６１０は、市販で入手可能なＯＳで動作する市販で入手可能な機械の演算部６１６を含む。ワークステーション６１０は、ＭＲＩ装置に走査指示の登録を可能とする操作インターフェイスを与える。ワークステーション６１０は、パルス系列サーバ６１８、データ取得サーバ６２０、データ処理サーバ６２２、及び、データ蓄積サーバ６２３の４つのサーバに接続されている。ワークステーション６１０及び各サーバ６１８、６２０、６２２、６２３は、互いに通信可能に接続されている。
【００４９】
パルス系列サーバ６１８は、ワークステーション６１０からダウンロードした命令に応じて作用し、傾斜装置６２４及びＲＦ装置６２６を操作する。所定の走査の実行に必要な傾斜波形が形成され、傾斜装置６２４に適用される。傾斜装置６２４は、装置６２８内の傾斜コイルを励起し、磁場傾斜Ｇｘ、Ｇｙ、Ｇｚを生成し、ＭＲ信号の位置エンコードに使用される。傾斜コイル装置６２８は、分極磁石６３２及びＲＦコイル６３４全体を含む磁石装置６３０の部分を形成する。
【００５０】
ＲＦ励起波形は、ＲＦ装置６２６によりＲＦコイル６３４に加えられ、所定の磁気共鳴パルス系列を実施する。ＲＦコイル６３４又は別個の局在コイル（図６に図示せず）で検出される応答ＭＲ信号は、パルス系列サーバ６１８で形成されるコマンドの指令で、ＲＦ装置６２６での受け取り、増幅、復調、フィルター及びデジタル化される。ＲＦ装置６２６は、ＭＲパルス系列に使用される幅広いＲＦパルスを形成するＲＦ送信器を含む。ＲＦ送信器は、パルス系列サーバ６１８からの走査指示及び走査指令に応じ、所望の周波数、位相及びパルス振幅波形を形成する。生成されたＲＦパルスは、ＲＦコイル全体、又は、１つ又はそれ以上の局在コイル又はコイル配列（図６に図示せず）に加えられる。
【００５１】
ＲＦ装置６２６は、１つ又はそれ以上のＲＦ受信路も含んでいる。各ＲＦ受信路は、接続されたコイルで受信されるＭＲ信号を増幅するＲＦ増幅器と、受信ＭＲ信号のＩ及びＱ直交成分を検出及びデジタル化する検出器とを含んでいる。受信されたＭＲ信号の大きさは、採取点でのＩ及びＱ成分の二乗の和の平方根、
Ｍ＝√(Ｉ２＋Ｑ２)
及び、受信ＭＲ信号の位相で決定される。
φ＝ｔａｎ−１(Ｑ/Ｉ)
【００５２】
パルス系列サーバ６１８は、生理学的取得制御部６３６から患者データを適宜受け取る。制御部６３６は、患者に接続された幾つかの異なるセンサ（例、電極からのＥＣＧ信号又はベローズからの呼吸信号）からの信号を受け取る。これらの信号は、パルス系列サーバ６１８で使用され、被検者の呼吸又は心臓の鼓動で走査の実行を同期又はゲートさせる。
【００５３】
パルス系列サーバ６１８は、走査室インターフェイス回路６３８にも接続している。インターフェイス回路６３８は、患者の状態及び磁気装置に関連する種々のセンサから信号を受け取る。走査室インターフェイス回路６３８を経て、患者位置決め装置６４０は、走査中、患者を所望の位置に動かすコマンドを受け取る。
【００５４】
ＲＦ装置６２６で形成されたデジタルＭＲ信号サンプルは、データ取得サーバ６２０で受信される。データ取得サーバ６２０は、ワークステーション６１０からダウンロードされる指示に呼応して、実時間ＭＲデータを受け取り、過剰データによるデータ損失を防ぐバッファ蓄積部を提供する。幾つかの走査で、データ取得サーバ６２０は、取得したＭＲデータをデータ処理サーバ６２２に送るだけの場合がある。しかし、取得したＭＲデータから誘導される情報で、走査の更なる実行を制御する走査の場合、データ取得サーバ６２０は、そのような情報を形成するようにプログラムされ、情報をパルス系列サーバ６１８に運ぶ。例えば、予備走査中、ＭＲデータが取得され、そのＭＲデータは、パルス系列サーバ６１８で実行されるパルス系列の較正に使用される。又、自動信号を走査中に取得させ、ＲＦ又は傾斜装置の操作変数の調整、又は、ｋ空間が採取される視野順の制御に使用できる。データ取得サーバ６２０は、磁気共鳴血管造影（ＭＲＡ）走査での造影剤の到達の検出に使用されるＭＲ信号の処理に使用しても良い。全ての実施例で、データ取得サーバ６２０は、ＭＲデータを取得し、そのデータを実時間で処理し、走査の制御に使用される情報を形成する。
【００５５】
データ処理サーバ６２２は、データ取得サーバ６２０からＭＲデータを受け取り、ワークステーション６１０からダウンロードされた指示に従い、そのデータを処理する。処理は、例えば、２又は３次元画像を形成するための未処理ｋ空間ＭＲデータのフーリエ変換、フィルタの再構成画像への適用、取得ＭＲデータの逆投影法画像再構成の実施、機能ＭＲ画像の計算、及び、動画又は流れる画像の計算等を含む。
【００５６】
データ処理サーバ６２２で再構成された画像は、ワークステーション６１０に戻されて蓄積される。実時間画像が、データベースメモリキャッシュ（図示せず）に蓄積される。キャッシュから、磁石装置６３０近傍に位置する操作表示部６１２又は表示部６４２に向けて、担当医の使用に出力される。バッチモード画像又は選択された実時間画像は、ディスク蓄積部６４４上のホストデータベースに蓄積される。画像が再構成され、蓄積部に転送されると、データ処理サーバ６２２は、ワークステーション６１０のデータ蓄積サーバ６２３に通知する。操作者は、画像の保管、フィルムの形成、又は、画像をネットワーク経由で他の施設に送るためにワークステーション６１０を使用する。
【００５７】
従来のＭＲＩ装置を本発明で使用できるが、１つの実施形態では、医者の接近を可能にする設計にしたＭＲＩ装置を使用する。特に図７に示すように、術中ＭＲ画像処理が実施される時、患者は、磁石装置６３０内に配置され、患者の関与領域を、２つの離れた磁石の輪６９０、６９２間に位置する装置の等角点近傍に位置させる。磁石の輪６９０、６９２間に立つ医師は、患者の重要な領域への接近に制限されない。処置中、患者は動いても良く、医師が動かしても良い。集束超音波（ＦＵＳ）切除装置６８０に接続された超音波振動子（図示せず）を医師が操作し、患者の選択した領域を治療しても良い。表示部６４２に形成される画像を、医師が、切除装置の支援及び治療終了の決定に使用できる。
【００５８】
画像再構成
【００５９】
図８を参照されたい。本発明の実施形態は、画像装置で画像データの取得で開始する（ステップ８００）。図４Ａ、４Ｂに示すようなＣアームＸ線画像装置で実施する場合、画像データは、Ｘ線投影空間を採取する一連のＸ線投影として取得される。図６、７に示すようなＭＲＩ装置で実施する場合、画像データは、ｋ空間を採取する一連のｋ空間投影として取得される。１つの実施形態として、取得された画像データは、時間分解画像データで、画像フレームの時系列に対応する。より具体的には、「画像フレーム」は、対応する時点で取得される画像データ又は「タイムフレーム」から再構成される。
【００６０】
取得された画像データから、次に、初期画像フレームが、各々の対応するタイムフレームに対し再構成される（ステップ８０２）。これは標準的な画像再構成である。例えば、画像データがＸ線画像データの場合、フィルタ補正逆投影の画像再構成法が使用される。
例えば、画像データが一連のｋ空間投影の場合、逆投影法も使用される。他の方法として、ｋ空間投影データを直交座標格子に再編し、従来のフーリエ変換再構成法の使用が出来る。各タイムフレームで取得される画像データは、下取りの採取で、画像の縞、画像のギザギザ等の画像の不自然さが、上記の再構成された初期画像フレームに存在する。これら初期画像フレームは、被検者の先行画像の形成に使用される（ステップ８０４）。例えば、先行画像は、入手できる初期画像フレームの平均化で形成される。
【００６１】
続いて、得られた先行画像、取得画像データ、及び、現在入手できる初期画像フレームを使用して、対応する画像フレームの推定値を形成する（ステップ８０６）。画像フレーム推定値を、例えば、幾何学的再構成技術（ＡＲＴ）を使用して形成する。これらの画像フレーム推定値を利用して、先行画像と共に、上記記載の最小化法で被検者の最終画像フレームを再構成する（ステップ８０８）。特に、次の目的関数が反復して最小化される。
||Ψ１(Ｉａｖｅ−Iｐ)||１＋α||Ψ２Ｉｔ||１＋βΣｔ＝２Ｎ||Ψ３(Iｔ−Iｔ−１)||１
Ψ１、Ψ２、Ψ３：スパース変換
Ｉａｖｅ：現在入手可能画像フレーム推定値Ｉｔ
Ｉｐ：先行画像
α、β：正則化変数
Ｉｔ：ｔｔｈの画像フレーム
上覧の目的関数は、式（２）、式（７）、式（１０）と少なくとも２点で異なる。第１に、目的関数の最後の項は、
βΣｔ＝２Ｎ||Ψ３(Ｉｔ−Ｉｔ−１)||１
入手可能な画像フレーム推定値の全てに及ぶ２つの隣接する画像フレーム推定値間の差異の和である。このように、現在の画像フレームＩｔに対する最近接画像フレームＩｔ−１は、
式（２）、式（７）、式（１０）の先行画像と同様に作用する。下取りデータから現在の画像フレームの再構成を制約するために使用される。第２として、目的関数の第１項は、
||Ψ１(Ｉａｖｅ−Iｐ)||１
先行画像と入手可能な画像フレーム推定値の平均との差異を計算する。この理由から、前述の目的関数を用いる本発明の実施形態は、実時間画像の適用に適している。実時間画像の場合、データは連続して取得され、データの取得中に画像が再構成される。この方法は、例えば、Ｘ線又はＭＲ画像装置で誘導される治療処置中で取得される画像の品質の向上を可能にする。
【００６２】
適切な停止基準に達すると（ステップ８１０）、画像再構成法は中止され、所望の画像フレームが蓄積又はモニタに表示される。停止基準が満たされない場合、画像フレーム推定値は、目的関数の最小化の出力で更新される（ステップ８１２）。続いて、これら更新された画像フレーム推定値は、上覧の目的関数で利用される。典型的な停止基準は、２つの隣接する繰返しの差異が前以て設定された閾値よりも小さい場合である。
【００６３】
本発明は、１つ又はそれ以上の実施形態で記載している。これら記載したものの等価物、変更、改変は可能で、本発明の範囲内になる。
【特許請求の範囲】
【請求項１】
被検者の時系列画像フレームを画像装置で再構成する方法が、
ａ）前記画像装置で前記被検者から画像データを取得するステップと、
ｂ）前記取得した画像データを使用して先行画像を形成するステップと、
ｃ）前記取得したデータと前記先行画像とを使用して複数の画像フレーム推定値を形成するステップと、
ｄ）前記先行画像と前記複数の画像フレーム推定値とを使用して複数の画像フレームを再構成するステップと、
ｅ）停止基準が満たされるまでステップｄ）を反復して繰り返すステップと、を含む画像再構成の方法において、
前記複数の画像フレーム推定値が、各反復後で更新されることを特徴とする画像再構成の方法。
【請求項２】
前記ステップｄ）が、現在のタイムフレームからの前記画像フレーム推定値と前のタイムフレームからの前記画像フレーム推定値との差異を計算することを含んでいることを特徴とする請求項１記載の画像再構成の方法。
【請求項３】
前記ステップｂ）が、複数の初期画像フレームを再構成し、前記複数の初期画像フレームの平均化で前記先行画像を形成することを特徴とする請求項１記載の画像再構成の方法。
【請求項４】
前記画像装置が、Ｘ線蛍光画像装置、ＣアームＸ線画像装置、磁気共鳴画像（ＭＲＩ）装置、及び、治療ＭＲＩ装置の少なくとも１つであることを特徴とする請求項１記載の画像再構成の方法。
【請求項５】
前記ステップｃ）が、代数的再構成技術（ＡＲＴ）を実行することを特徴とする請求項１記載の画像再構成の方法。
【請求項６】
前記ステップｄ）が、
||Ψ１(Ｉａｖｅ−Iｐ)||１＋α||Ψ２Ｉｔ||１＋βΣｔ＝２Ｎ||Ψ３(Iｔ−Iｔ−１)||１
の目的関数を最小化するステップを含むことを特徴とする請求項１記載の画像再構成の方法。
Ψ１：第１スパース変換
Ψ２：第２スパース変換
Ψ３：第３スパース変換
Ｉａｖｅ：前記複数の画像フレーム推定値の平均
Ｉｐ：前記先行画像
α：第１正則化変数
Ｉｔ：現在のタイムフレームに対する画像フレーム推定値
β：第２正則化変数
Ｉｔ−１：前のタイムフレームに対する画像フレーム推定値
【請求項７】
画像装置で被検者の画像を再構成する方法が、
ａ）前記画像装置で前記被検者から画像データを取得するステップと、
ｂ）前記取得した画像データを使用して前記被検者の先行画像を形成するステップと、
ｃ）前記取得したデータと前記先行画像とを使用して現在のタイムフレームに対応する現在の画像フレーム推定値を形成するステップと、
ｄ）新しいタイムフレームを選択し、前記現在の画像フレーム推定値を前の画像フレーム推定値として蓄積するステップと、
ｅ）前記新しいタイムフレームに対応する前記現在の画像フレーム推定値を形成するステップｃ）を繰り返すステップと、
ｆ）前記前の画像フレーム推定値、前記新しい現在の画像フレーム推定値、及び、前記先行画像を使用して、前記新しいタイムフレームに対応する画像フレームを再構成するステップとを備えていることを特徴とする画像再構成の方法。
【請求項８】
前記ステップｄ）〜ｆ）が、複数のタイムフレームに対応する複数の画像フレームを形成するように繰り返されていることを特徴とする請求項７記載の画像再構成の方法。
【請求項９】
前記ステップｂ）が、複数の初期画像フレームを再構成し、前記複数の初期画像フレームを平均化して前記先行画像を形成することを特徴とする請求項７記載の画像再構成の方法。
【請求項１０】
前記ステップｃ）が、代数的再構成技術（ＡＲＴ）の実行を含んでいることを特徴とする請求項記７載の画像再構成の方法。
【請求項１１】
前記画像装置が、Ｘ線蛍光画像装置、ＣアームＸ線画像装置、磁気共鳴画像（ＭＲＩ）装置及び治療ＭＲＩ装置の少なくとも１つであることを特徴とする請求項７記載の画像再構成の方法。
【請求項１】
被検者の時系列画像フレームを画像装置で再構成する方法が、
ａ）前記画像装置で前記被検者から画像データを取得するステップと、
ｂ）前記取得した画像データを使用して先行画像を形成するステップと、
ｃ）前記取得したデータと前記先行画像とを使用して複数の画像フレーム推定値を形成するステップと、
ｄ）前記先行画像と前記複数の画像フレーム推定値とを使用して複数の画像フレームを再構成するステップと、
ｅ）停止基準が満たされるまでステップｄ）を反復して繰り返すステップと、を含む画像再構成の方法において、
前記複数の画像フレーム推定値が、各反復後で更新されることを特徴とする画像再構成の方法。
【請求項２】
前記ステップｄ）が、現在のタイムフレームからの前記画像フレーム推定値と前のタイムフレームからの前記画像フレーム推定値との差異を計算することを含んでいることを特徴とする請求項１記載の画像再構成の方法。
【請求項３】
前記ステップｂ）が、複数の初期画像フレームを再構成し、前記複数の初期画像フレームの平均化で前記先行画像を形成することを特徴とする請求項１記載の画像再構成の方法。
【請求項４】
前記画像装置が、Ｘ線蛍光画像装置、ＣアームＸ線画像装置、磁気共鳴画像（ＭＲＩ）装置、及び、治療ＭＲＩ装置の少なくとも１つであることを特徴とする請求項１記載の画像再構成の方法。
【請求項５】
前記ステップｃ）が、代数的再構成技術（ＡＲＴ）を実行することを特徴とする請求項１記載の画像再構成の方法。
【請求項６】
前記ステップｄ）が、
||Ψ１(Ｉａｖｅ−Iｐ)||１＋α||Ψ２Ｉｔ||１＋βΣｔ＝２Ｎ||Ψ３(Iｔ−Iｔ−１)||１
の目的関数を最小化するステップを含むことを特徴とする請求項１記載の画像再構成の方法。
Ψ１：第１スパース変換
Ψ２：第２スパース変換
Ψ３：第３スパース変換
Ｉａｖｅ：前記複数の画像フレーム推定値の平均
Ｉｐ：前記先行画像
α：第１正則化変数
Ｉｔ：現在のタイムフレームに対する画像フレーム推定値
β：第２正則化変数
Ｉｔ−１：前のタイムフレームに対する画像フレーム推定値
【請求項７】
画像装置で被検者の画像を再構成する方法が、
ａ）前記画像装置で前記被検者から画像データを取得するステップと、
ｂ）前記取得した画像データを使用して前記被検者の先行画像を形成するステップと、
ｃ）前記取得したデータと前記先行画像とを使用して現在のタイムフレームに対応する現在の画像フレーム推定値を形成するステップと、
ｄ）新しいタイムフレームを選択し、前記現在の画像フレーム推定値を前の画像フレーム推定値として蓄積するステップと、
ｅ）前記新しいタイムフレームに対応する前記現在の画像フレーム推定値を形成するステップｃ）を繰り返すステップと、
ｆ）前記前の画像フレーム推定値、前記新しい現在の画像フレーム推定値、及び、前記先行画像を使用して、前記新しいタイムフレームに対応する画像フレームを再構成するステップとを備えていることを特徴とする画像再構成の方法。
【請求項８】
前記ステップｄ）〜ｆ）が、複数のタイムフレームに対応する複数の画像フレームを形成するように繰り返されていることを特徴とする請求項７記載の画像再構成の方法。
【請求項９】
前記ステップｂ）が、複数の初期画像フレームを再構成し、前記複数の初期画像フレームを平均化して前記先行画像を形成することを特徴とする請求項７記載の画像再構成の方法。
【請求項１０】
前記ステップｃ）が、代数的再構成技術（ＡＲＴ）の実行を含んでいることを特徴とする請求項記７載の画像再構成の方法。
【請求項１１】
前記画像装置が、Ｘ線蛍光画像装置、ＣアームＸ線画像装置、磁気共鳴画像（ＭＲＩ）装置及び治療ＭＲＩ装置の少なくとも１つであることを特徴とする請求項７記載の画像再構成の方法。
【図１】
【図２】
【図３】
【図４】
【図５Ａ】
【図５Ｂ】
【図６】
【図７】
【図８】
【図２】
【図３】
【図４】
【図５Ａ】
【図５Ｂ】
【図６】
【図７】
【図８】
【公表番号】特表２０１１−５０９７１７（Ｐ２０１１−５０９７１７Ａ）
【公表日】平成２３年３月３１日（２０１１．３．３１）
【国際特許分類】
【出願番号】特願２０１０−５４２４３１（Ｐ２０１０−５４２４３１）
【出願日】平成２１年１月１４日（２００９．１．１４）
【国際出願番号】ＰＣＴ／ＵＳ２００９／０３０９９４
【国際公開番号】ＷＯ２００９／０９１８２４
【国際公開日】平成２１年７月２３日（２００９．７．２３）
【出願人】（３９００２３６４１）ウイスコンシン アラムナイ リサーチ ファウンデーシヨン (61)
【氏名又は名称原語表記】ＷＩＳＣＯＮＳＩＮ ＡＬＵＭＮＩ ＲＥＳＥＡＲＣＨ ＦＯＵＮＤＡＴＩＯＮ
【Ｆターム（参考）】
【公表日】平成２３年３月３１日（２０１１．３．３１）
【国際特許分類】
【出願日】平成２１年１月１４日（２００９．１．１４）
【国際出願番号】ＰＣＴ／ＵＳ２００９／０３０９９４
【国際公開番号】ＷＯ２００９／０９１８２４
【国際公開日】平成２１年７月２３日（２００９．７．２３）
【出願人】（３９００２３６４１）ウイスコンシン アラムナイ リサーチ ファウンデーシヨン (61)
【氏名又は名称原語表記】ＷＩＳＣＯＮＳＩＮ ＡＬＵＭＮＩ ＲＥＳＥＡＲＣＨ ＦＯＵＮＤＡＴＩＯＮ
【Ｆターム（参考）】
[ Back to top ]