血流動態解析装置、血流動態解析プログラム、流体解析装置及び流体解析プログラム
【課題】処理時間を増大させることなく流体の動態に関する情報を安定に算出することが可能な血流動態解析装置を提供することである。
【解決手段】 本発明の実施形態に係る血流動態解析装置は、組織及び動脈の時間濃度曲線算出手段、第１及び第２の構成手段、逆畳み込み積分手段並びに血流情報算出手段を備える。組織の時間濃度曲線算出手段は、造影画像データの時相間における画素値変化を計測して組織の時間濃度曲線を算出する。動脈の時間濃度曲線算出手段は、組織の時間濃度曲線から動脈の時間濃度曲線を算出する。第１の構成手段は、組織の時間濃度曲線から第１の集合を構成する。第２の構成手段は、動脈の時間濃度曲線から第２の集合を構成する。逆畳み込み積分手段は、第１及び第２の集合の逆畳み込み積分を行い、組織の伝達関数を算出する。血流情報算出手段は、伝達関数に基づいて、血流動態に関する情報を算出する。
【解決手段】 本発明の実施形態に係る血流動態解析装置は、組織及び動脈の時間濃度曲線算出手段、第１及び第２の構成手段、逆畳み込み積分手段並びに血流情報算出手段を備える。組織の時間濃度曲線算出手段は、造影画像データの時相間における画素値変化を計測して組織の時間濃度曲線を算出する。動脈の時間濃度曲線算出手段は、組織の時間濃度曲線から動脈の時間濃度曲線を算出する。第１の構成手段は、組織の時間濃度曲線から第１の集合を構成する。第２の構成手段は、動脈の時間濃度曲線から第２の集合を構成する。逆畳み込み積分手段は、第１及び第２の集合の逆畳み込み積分を行い、組織の伝達関数を算出する。血流情報算出手段は、伝達関数に基づいて、血流動態に関する情報を算出する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明の実施形態は、血流動態解析装置、血流動態解析プログラム、流体解析装置及び流体解析プログラムに関する。
【背景技術】
【０００２】
従来、Ｘ線CT（Computed Tomography）装置やMRI（Magnetic Resonance Imaging）装置では、静脈から造影剤を投与し、時系列画像データを収集して、その画像を元に解析し、組織内の血流動態に関する情報を得る検査が行われている。これはパーヒュージョン検査と呼ばれ、造影剤の撮影断面への集中の度合いが画像の濃度変化として収集できることを利用したものである。
【０００３】
例えば、肺循環や造影剤投与のばらつきをなくすために、組織に流入する動脈の時間濃度曲線（Time Concentration Curve：TCC）：ａｉを入力関数として、測定した組織の時間濃度曲線Ｃｉとのデコンボリューション（Deconvolution：逆畳み込み積分）を行い、組織固有のインパルス応答関数（伝達関数）Ｒｉ（residue function）を求め、得られたインパルス応答関数から、血流動態を定量的に表す指標である血流量（脳の場合、CBF：Cerebral Blood Flow）、平均通過時間：MTT（Mean Transit Time）、および、血液量（脳の場合、CBV：Cerebral Blood Volume）などのパラメータを算出している。この解析方法は、デコンボリューション法と呼ばれている。
【０００４】
尚、時間濃度曲線は、造影剤濃度に比例する量の時間変化を表す曲線である。Ｘ線CT装置で収集された造影画像データから求められる造影剤濃度に比例する量の時間変化曲線は、TDC (Time Density Curve)とも呼ばれるが、ここではTCCと表記する。
【０００５】
例えば、非特許文献１では、デコンボリューションに、standard SVD（Singular Value Decomposition）法やblock circulant SVD法を適用している。
【０００６】
ここで、standard SVD法について説明する。まず、収集された画像データにおいて、画素ｊ（位置ｒｊ）の画素値の変化を計測して組織の時間濃度曲線Ｃｉ（ｒｊ）を算出する。また、算出した組織の時間濃度曲線Ｃｉから、動脈の時間濃度曲線Ａｉを次式（１）に示すようにして算出する。式（１）において、ｉは、フェーズ番号を表わし、ΩＡは、動脈に相当する画素の集合を表わしている。
【数１】
【０００７】
組織の時間濃度曲線Ｃｉが、動脈の時間濃度曲線Ａｉとインパルス応答関数Ｘのデコンボリューション（畳み込み積分）であるとすると、次式（２）に示すようなシステム方程式で表わすことができる。式（２）において、ＴＲは、画像撮影の時間間隔（繰り返し時間 TR: repetition time)を表わし、Ｒｊは、インパルス応答関数を表わしている。
【数２】
【０００８】
例えば、動脈の時間濃度曲線ａｉおよび組織の時間濃度曲線Ｃｉを時間軸方向に離散化したもの、つまり、ｉｂｓ≦ｉ≦Ｌ−１のフェーズ範囲をＬ個の等区間に分割し、Ｌ個のフェーズのデータを解析に用いるとすると、上記の式（２）は、行列形式でＣ＝ＡＸとなり、次式（３）に示すように表わすことができる。式（３）において、Ｃは、Ｌ×１の列ベクトルであり、Ａは、Ｌ×Ｌの係数行列であり、Ｘは、Ｌ×１の列ベクトルである。
【数３】
【０００９】
上記の式（３）において、ｉｂｓは、解析対象となるフェーズ範囲の先頭フェーズを表わしており、動脈および組織の造影濃度が上昇し始めるフェーズ、もしくは、それよりも前のフェーズを用いる。また、係数行列の各要素値ａi,jは、次式（４）により表わされる。式（４）において、ｉ＝０，１，・・・Ｌ−１であり、ｊ＝０，１，...Ｌ−１である。
【数４】
そして、上記の式（３）の近似解を、特異値分解により求める。すなわち、行列Ａを次式（５）に示すように３つの行列に分解することができる。
［数５］
Ａ＝ＵＤＶＴ ・・・（５）
【００１０】
上記の式（５）において、Ｕは、ｍ行ｋ列の正規直交行列であり、Ｄは、ｋ行ｋ列の対角行列であり、Ｖは、ｎ行ｋ列の正規直交行列である。上記の式（３）に用いる場合、ｍ＝ｎ＝ｋ＝Ｌである。対角行列とは、対角成分以外の成分がすべて０である行列であり、正規直交行列とは、対角成分が１である対角行列である。
【００１１】
ここで、Ｕ＝（ｕ１，ｕ２，...ｕｋ）、Ｖ＝（ｖ１，ｖ２，...ｖｋ）、Ｄ＝diag（ｗ０，ｗ１，...ｗｋ）であるとき、ｗｋを降順に並び替え、ｕｋ、ｖｋもその順序に対応するように入れ替えを行っておく。行列Ａの特異値分解の結果を用いて、Ｘの最小２乗解である最小ノルム推定値を次式（６）により算出することができる。
【数６】
【００１２】
上記の式（６）において、Ｗは、Ｄの対角要素（特異値）の逆数からなる対角行列であって、Ｗ＝diag（１／ｗｏ，１／ｗ１，...１／ｗＧ−１，０，０，...）である。ここで、最大特異値のＰｓｖｄ倍に満たない特異値に対しては、逆数ではなく０を対角要素とする。従って、Ｇ個の対角要素だけが０でない値を持つようになる。Ｐｓｖdは、正則化の強さを左右するパラメータであり、大きい値ほど強い正則化が行われる。
そして、ノルム推定値が求められたら、次式（７）に示すようにしてインパルス応答関数Ｒｉを算出することができる。
【数７】
上記の式（７）で算出されたインパルス応答関数Ｒｉから、組織血流量、組織血液量、および平均通過時間を算出することができる。
【００１３】
ところで、血流は、組織の毛細血管へ流れるが、造影剤が組織に到達する時間は、組織の場所によって異なり、設定した動脈領域から求めた動脈の時間濃度曲線と各画素の組織の時間濃度曲線は、厳密には時間軸が一致しない。従って、極端な場合、動脈の造影濃度が上昇する前に（動脈の時間濃度曲線が立ち上がる前に）、組織の造影濃度が上昇してしまう（組織の時間濃度曲線が先に立ち上がってしまう）こととなる。
【００１４】
上述したstandard SVD法では、このような時間方向のずれ（ディレイ）がないことを仮定としている。しかしながら、動脈の時間濃度曲線と組織の時間濃度曲線の間には、必ず時間方向のずれが存在するため、ディレイ時間の値によってデコンボリューション結果が変わり、算出される組織血流量の値が影響を受けてしまうこととなる。
【００１５】
そこで、block circulant SVD法では、時間方向のずれに対する依存性を軽減するために、次式（８）に示すようなブロック循環行列であるシステム行列を用いるようにする。
【数８】
【００１６】
このblock circulant SVD法において、基本的な処理は、standard SVD法と同様であり、行列形式Ｃ＝ＡＸであるが、Ａが、２Ｌ×２Ｌの係数行列となっている。そのため、動脈の時間濃度曲線を構成する係数行列Ａの各要素値の後方に０が追加され、インパルス応答関数を構成する列ベクトルＸの各要素値の後方に０が追加され、組織の時間濃度曲線を構成する列ベクトルＣの各要素値の後方に０が追加され、使用する測定データの２倍の要素数となるように拡張される。特に、係数行列Ａは、行および列がいずれも２倍、もしくはそれ以上に拡張され、さらに列ベクトルの要素が周期的になるように配列し直されている。
【００１７】
このように、システム行列を上記の式（８）に示すように構成することで、時間方向のずれが存在した場合、インパルス応答関数が時間軸方向にシフトし、その形状が相似形を保つようになる。従って、インパルス応答関数を時間軸方向にシフトさせることによって、時間方向のずれに対する依存性を小さくした、組織血流量、組織血液量、および平均通過時間を算出することができる。
【００１８】
しかしながら、システム行列において、動脈の時間濃度曲線を構成する係数行列Ａが２×２倍、組織の時間濃度曲線を構成する列ベクトルが２倍となることから、解析処理時間は８倍以上にも増大してしまい、実用的ではない。
【００１９】
そこで、非特許文献２では、解析処理時間を増大させずに、standard SVD法のシステム行列を応用している。すなわち、非特許文献２で提案されている方法では、システム行列として上記の式（４）を用い、既知ベクトルには、先頭フェーズからオフセット位置までの組織の時間濃度曲線の要素をノイズ（信号のひずみ）として計算上から除外している。ここでは、上記の式（３）において、ベクトルＣからオフセット位置までの先頭数フェーズが除かれることから、後方には、同数の０の要素が追加されることになる。
このように、非特許文献２の技術では、システム行列を拡張せずに、時間方向のずれに対する依存性を小さくするようにしている。
【先行技術文献】
【非特許文献】
【００２０】
【非特許文献１】Ona Wu et al,"Tracer Arrival Timing-Insensitive Technique for Estimating Flow in MR Perfusion-Weighted Imaging Using Singular Value Decomposition With a Block-Circulant Deconvolution Matrix",Magnetic Resonance in Medicine 50:164-174(2003)
【非特許文献２】M.R.Smith, H.Lu, S.Trochet, and R.Frayne "Removing the Effect of SVD Algorithmic Artifacts Present in Quantitative MR Perfusion Studies", Magnetic Resonance in Medicine 51:631-634(2004)
【発明の概要】
【発明が解決しようとする課題】
【００２１】
非特許文献１の技術において、standard SVD法では、時間方向のずれを考慮していないため、デコンボリューション結果の正確性に欠ける。一方、block circulant SVD法では、時間方向のずれを考慮しているものの、解析処理時間が膨大であり、いずれも実用的ではない課題があった。
【００２２】
特に、Ｘ線CT装置やMRI装置の運用では、通常、多数の患者の検査を実施しなければならない。また、画像撮影後、直ちに解析処理を行って、技師が、処理結果を見て撮影に不備があるか否かを確認しなければならない。従って、解析処理時間が膨大になると、実施可能な検査数が減少してしまうとともに、患者にとっては待ち時間が増え、病院にとっては患者獲得の減少になり、双方に不利益をもたらすこととなる。
【００２３】
そこで、非特許文献２の技術では、解析処理時間を増大させずに処理を行っているが、先頭フェーズからオフセット位置までの組織の時間濃度曲線の要素をノイズ（信号のひずみ）として計算上から除外してしまうことから、オフセット位置に結果が依存してしまうという新たな問題が生じてしまう。しかしながら、オフセット位置を正しく同定することは難しいため、結果的に血流動態に関する情報を安定に求めることができないという課題があった。
【００２４】
本発明はこのような状況に鑑みてなされたものであり、その目的は、処理時間を増大させることなく流体の動態に関する情報を安定に算出することが可能な血流動態解析装置、血流動態解析プログラム、流体解析装置及び流体解析プログラムを提供することである。
【課題を解決するための手段】
【００２５】
本発明の実施形態に係る血流動態解析装置は、組織の時間濃度曲線算出手段、動脈の時間濃度曲線算出手段、第１の構成手段、第２の構成手段、逆畳み込み積分手段及び血流情報算出手段を備える。組織の時間濃度曲線算出手段は、医用画像撮影装置により得られた複数の時相における被検体の造影画像データに基づいて、前記複数の時相間における各画素の画素値の変化を計測して組織における造影剤の濃度に対応する量の時間変化を表す時間濃度曲線を算出する。動脈の時間濃度曲線算出手段は、前記組織に対応する前記時間濃度曲線から動脈領域を設定し、動脈に対応する時間濃度曲線を算出する。第１の構成手段は、前記組織に対応する前記時間濃度曲線を、第１の時刻を基準にして離散化し、異なる時刻での前記組織における前記造影剤の濃度に対応する値を含む複数の要素を１次元に配列した第１の集合を構成する。第２の構成手段は、前記動脈に対応する前記時間濃度曲線を、前記第１の時刻より後方の第２の時刻を基準にして離散化し、異なる時刻での前記動脈における前記造影剤の濃度に対応する値を含む複数の要素を２次元に配列した第２の集合を構成する。逆畳み込み積分手段は、前記第１の集合と前記第２の集合の逆畳み込み積分を行い、前記組織の伝達関数を算出する。血流情報算出手段は、前記伝達関数に基づいて、血流動態に関する情報を算出する。
【００２６】
また、本発明の実施形態に係る血流動態解析プログラムは、コンピュータを、上述の組織の時間濃度曲線算出手段、動脈の時間濃度曲線算出手段、第１の構成手段、第２の構成手段、逆畳み込み積分手段及び血流情報算出手段として機能させる。
【００２７】
また、本発明の実施形態に係る流体解析装置は、濃度変化算出手段及び流体情報算出手段を備える。濃度変化算出手段は、被検体の複数のフレームの画像データに基づいて、組織に対応する画素において画素値の濃度に応じた量の前記複数のフレーム間における第１の変化を算出する一方、流体に対応する画素値の濃度に応じた量の前記複数のフレーム間における第２の変化を算出する。流体情報算出手段は、前記第１の変化に基づいて前記第１の変化の第１の点に対応する要素を最初の要素とする解析対象となる数の第１の複数の要素を作成する一方、前記第２の変化に基づいて前記第２の変化の前記第１の点よりも後の第２の点に対応する要素を最初の要素とする前記解析対象となる数の第２の複数の要素を作成し、前記第１の複数の要素及び前記第２の複数の要素に基づいて逆畳み込み積分処理を伴って流体の動態に関する情報を算出する。
また、本発明の実施形態に係る流体解析プログラムは、コンピュータを、上述の濃度変化算出手段及び流体情報算出手段として機能させる。
【図面の簡単な説明】
【００２８】
【図１】本発明の第１の実施形態に係る血流動態解析装置の構成例を示す図。
【図２】局所血流解析アプリケーションを実行するCPUの機能の一例を示すブロック図。
【図３】ピークフェーズを説明する図。
【図４】組織の血流動態に関する情報の算出処理を説明するフローチャート。
【図５】第１の実施形態に係るシステム行列を模式的に示す図。
【図６】本発明の第２の実施形態に係る流体解析装置の構成例を示す機能ブロック図。
【図７】図６に示すMRI装置のイメージング部において実行される非造影流体イメージング用の撮像条件の第１の例を示す図。
【図８】図６に示すMRI装置のイメージング部において実行される非造影流体イメージング用の撮像条件の第２の例を示す図。
【図９】図６に示すMRI装置のイメージング部において実行される非造影流体イメージング用の撮像条件の第３の例を示す図。
【図１０】図６に示すシステム行列作成部３ｄにおける解析対象となる要素の作成方法を説明する図。
【図１１】図６に示す流体解析装置におけるデータ処理の流れを示すフローチャート。
【発明を実施するための形態】
【００２９】
以下、本発明の実施の形態について図面を参照して詳細に説明する。
【００３０】
（第１の実施形態）
図１は、本発明の第１の実施形態に係る血流動態解析装置１の構成例を示す図である。
【００３１】
血流動態解析装置１は、CPU(Central Processing Unit)１ａ、ROM（Read Only Memory）１ｂ、RAM（Random Access Memory）１ｃ、および入出力インターフェイス１ｄが、バス１ｅを介して接続されている。入出力インターフェイス１ｄには、記憶部１ｆ、通信部１ｇ、表示部１ｈ、入力部１ｉ、およびリムーバブルメディア１ｊが接続されている。
【００３２】
CPU１ａは、入力部１ｉからの入力信号に基づいて血流動態解析装置１を起動するためのブートプログラムをROM１ｂから読み出して実行し、記憶部１ｆに格納されている各種オペレーティングシステムを読み出す。またCPU１ａは、入力部１ｉからの入力信号に基づいて各種の制御を行ったり、ROM１ｂや記憶部１ｆに記憶されたプログラムおよびデータを読み出してRAM１ｃにロードしたり、あるいはRAM１ｃから読み出されたプログラムのコマンドに基づいて、データ演算または加工などの一連の処理を実行する。
【００３３】
記憶部１ｆは、半導体メモリや磁気ディスクなどで構成されており、CPU１ａで実行されるプログラムやデータを記憶する。記憶部１ｆには、CPU１ａが実行するプログラムとして、例えば、局所血流動態に関する情報を解析するためのアプリケーションが用意される。以下、局所血流動態に関する情報を解析するこのアプリケーションを局所血流解析アプリケーションという。
【００３４】
通信部１ｇは、LAN（Local Area Network）カードやモデムなどで構成されており、血流動態解析装置１をLANやインターネットといった通信媒体に接続することを可能にする。すなわち通信部１ｇは、通信媒体から受信したデータを、入出力インターフェイス１ｄおよびバス１ｅを介してCPU１ａに送信し、CPU１ａからバス１ｅおよび入出力インターフェイス１ｄを介して受信したデータを、通信媒体に送信する。
【００３５】
表示部１ｈは、例えば液晶ディスプレイである。CPU１ａからバス１ｅおよび入出力インターフェイス１ｄを介して受信した信号に基づいて、CPU１ａの処理結果などを表示する。
【００３６】
入力部１ｉは、血流動態解析装置１の操作者が各種の操作を入力するキーボードやマウスなどの入力デバイスにより構成されている。操作者の操作に基づいて入力信号を生成し、入出力インターフェイス１ｄおよびバス１ｅを介してCPU１ａに送信する。
【００３７】
リムーバブルメディア１ｊは、例えば光ディスクやフレキシブルディスクなどである。図示せぬディスクドライブによって読み出されたデータが、入出力インターフェイス１ｄおよびバス１ｅを介してCPU１ａに送信され、CPU１ａからバス１ｅおよび入出力インターフェイス１ｄを介して受信したデータが、ディスクドライブによって書き込まれる。
【００３８】
第１の実施形態における血流動態解析装置１は、Ｘ線CT装置やMRI装置などの医用画像撮影装置（いわゆる医用モダリティ）により、ダイナミックス撮影により収集された画像データから血流動態を解析する装置である。このため、血流動態解析装置１は、かかる画像データを入手できる環境にあればよく、医用画像撮影装置と一体に構成されていてもよいし、医用画像撮影装置とは別体で構成されていてもよい。別体で構成される場合には、収集された画像データはリムーバブルメディア１ｊで、または通信部１ｇを介して図示せぬ医用画像撮影装置から血流動態解析装置１に送られる。
【００３９】
また、血流動態解析装置１は、血流動態解析結果を、通信部１ｇを介して、PACS(Picture Archiving and Communication System)に転送し、そこに保存させることができる。これによって、特定の画像データの血流動態解析結果の検索や閲覧などを容易に行うことが可能となる。
【００４０】
医用画像撮影装置には、Ｘ線CT装置やMRI装置のほか、超音波診断装置やPET（Positron Emission Tomography）装置などが含まれる。
図２は、局所血流解析アプリケーションを実行するCPU１ａの機能の一例を示すブロック図である。
【００４１】
図２に示すように、局所血流解析アプリケーションを実行するCPU１ａは、時系列データ取得部２ａ、造影濃度曲線算出部２ｂ、動脈の時間濃度曲線算出部２ｃ、システム行列作成部２ｄ、およびデコンボリューション部２ｅとして機能する。
【００４２】
時系列データ取得部２ａは、造影剤を静脈注入（ボーラス注入）した被検者について、医用画像撮影装置が所定時間撮影して生成した時系列の撮影データ（多時相の画像データ）を取得し、取得した時系列の撮影データを造影濃度曲線算出部２ｂに出力する。
造影濃度曲線算出部２ｂは、時系列データ取得部２ａから出力される時系列の撮影データを入力し、各画素の画素値の変化を計測して造影濃度曲線を算出する。
【００４３】
つまり、被検者に造影剤を比較的短時間で静脈注入すると、造影剤は静脈を通り心臓に戻り、肺を循環して再び心臓に戻り、動脈へ流入し、各臓器の組織に到達し、再び静脈を通り心臓へ戻る。ここで、例えば、医用画像撮影装置が脳の領域のCT画像やMRI画像を時系列で撮影していた場合、造影剤濃度の変化に応じた各画素の画素値の変化を計測することができる。造影濃度曲線算出部２ｂは、この計測した各画素の画素値の変化に基づいて、造影濃度曲線を算出することができる。
【００４４】
例えば、CT画像の場合、次式（９）に示すように、各時相の画素値Ｓｉから造影剤到着前の画素値Ｓ０を引き算することにより、造影剤濃度に比例する量（TDC：Time Density Curve）が算出される。式（９）において、ｉは、撮影時相を表している。
［数９］
Ｃｉ＝Ｓｉ−Ｓ０ ・・・（９）
【００４５】
また例えば、MRI画像の場合、次式（１０）に示すように、各時相の画素値Ｓｉと造影剤到着前の画素値Ｓ０の比の対数（ΔR2＊）から、造影剤濃度に比例する量（TCC：Time Concentration Curve）が算出される。式（１０）において、ｉは、撮影時相を表し、TEは、エコー時間を表している。
［数１０］
ΔＲ２＊＝−log（Ｓｉ／Ｓ０）／TE ・・・（１０）
【００４６】
造影濃度曲線算出部２ｂは、上記の式（９）または式（１０）で算出した造影濃度曲線を、組織の時間濃度曲線Ｃｉとし、算出値を動脈の時間濃度曲線算出部２ｃ、および、システム行列作成部２ｄに出力する。
【００４７】
この造影濃度曲線算出部２ｂは、複数の時相における被検体の造影画像データに基づいて、複数の時相間における各画素の画素値の変化を計測して組織における造影剤の濃度に対応する量の時間変化を表す時間濃度曲線を算出する組織の時間濃度曲線算出手段として機能する。
【００４８】
動脈の時間濃度曲線算出部２ｃは、造影濃度曲線算出部２ｂから出力される組織の時間濃度曲線Ｃｉを入力し、その曲線から、CT画像やMRI画像の解析単位となるボクセルでの動脈の時間濃度曲線を算出する。つまり、CT画像やMRI画像の解析単位となるボクセルの大きさは、０．５ｍｍ乃至３ｍｍ程度の大きさであり、組織の中の主要な動脈であれば、血管にほぼ包含されるボクセルが存在する。
【００４９】
動脈の時間濃度曲線算出部２ｃは、このようなボクセルを自動的あるいはユーザからの指定によって特定することで、動脈領域を設定し、その動脈領域での時間濃度曲線（ここでは、ｉ番目の撮影時刻での濃度値ａi）を算出する動脈の時間濃度曲線算出手段として機能する。なお、複数の画素が動脈領域として特定された場合には、例えば、その平均値をａiとする。
【００５０】
或いは、動脈の時間濃度曲線算出部２ｃは、必要に応じて、次式（１１）によるガンマ確率密度関数の定数倍の関数（ガンマ確率密度関数に比例する関数）を動脈の時間濃度曲線のカーブモデル（以下、「動脈カーブモデル」と称する）として適用し、算出した動脈の時間濃度曲線ａｉのカーブフィッティング（動脈カーブモデルへの当てはめ）を行う。通常、MRI画像ではカーブフィッティングを行う場合が多いが、CT画像ではカーブフィッティングを行わない場合が多い。
［数１１］
φ(t; t0, α, β, K)=K(t-t0)αexp{-β-1(t-t0)} ・・・（１１）
【００５１】
ここで言うガンマ確率密度関数とは、時間βに１回の割合で通過できる何重かの障壁αの全部を通過する所要時間の確率密度を表わしている。動脈関数は、肺の障壁を通過して分散された関数と考えることができるため、このガンマ確率密度関数を、動脈カーブモデルとして適用し、動脈カーブモデルへの当てはめ（近似）を行うようにする。
【００５２】
動脈の時間濃度曲線ａ（t）は、上記の式（１１）で表されるガンマ確率密度関数の定数倍の関数でのカーブフィッティングによって、次式（１２）により表わされる。式（１２）において、ｔは、時刻を表し、t0，α，β，Ｋは、いずれもパラメータを表している。αは、障壁の数であり、βは、障壁αの全部を通過する所要時間であり、Ｋは、線形パラメータである。
［数１２］
a(t)=φ(t; t0, α, β, K) ・・・（１２）
【００５３】
動脈の時間濃度曲線算出部２ｃは、上記の式（１２）における４つのパラメータt0，α，β，Ｋの全てを変数とし、次式（１３）に示す、動脈カーブモデルと動脈の濃度値ａｉとの２乗誤差を最小にするt0，α，β，Ｋを算出する。実際には、t0には、造影剤が対象臓器に到達する前の特定の時刻を用い、α，β，Ｋを算出する方法を用いる場合が多い。
【数１３】
動脈の時間濃度曲線算出部２ｃは、カーブフィッティングした動脈の時間濃度曲線ａｉの算出値を、システム行列作成部２ｄに出力する。
【００５４】
システム行列作成部２ｄは、造影濃度曲線算出部２ｂから出力される組織の時間濃度曲線Ｃｉの値を入力するとともに、動脈の時間濃度曲線算出部２ｃから出力される動脈の時間濃度曲線ａｉの値を入力する。
【００５５】
システム行列作成部２ｄは、入力した組織の時間濃度曲線Ｃｉおよび動脈の時間濃度曲線ａｉから、解析対象となる先頭フェーズｉｂｓおよびフェーズ範囲ｉｂｓ≦ｉ≦Ｌ−１を決定する。解析対象となる先頭フェーズｉｂｓは、動脈および組織の造影濃度が上昇し始めるフェーズ、もしくは、それよりも前のフェーズを用いる。フェーズとは、解析上のサンプリング間隔により決まる時刻のことであり、フェーズ範囲とはいわゆる解析範囲のことである。
【００５６】
システム行列作成部２ｄは、決定したフェーズ範囲における動脈の時間濃度曲線ａｉを時間軸方向に離散化したもの、つまり、ｉｂｓ≦ｉ≦Ｌ−１のフェーズ範囲をＬ個の等区間に分割したＬ個のデータを、次式（１４）に示すようにして算出する。
［数１４］
Ａｉ＝ａ（ｉＴＲ） ・・・（１４）
【００５７】
上記の式（１４）において、ＴＲは、画像撮影の時間間隔（繰り返し時間）を表わしているが、これに限らず、ダイナミック時相、または、ダイナミック時相から得られた実際の解析時相でもよい。また、測定値ａｉをそのままＡｉとして使用するようにしてもよい。
【００５８】
システム行列作成部２ｄは、組織の時間濃度曲線Ｃｉを時間軸方向に離散化することにより得た、異なる時刻での組織の濃度値からなる集合の、時間軸方向の前方に、Ｈ個の要素を追加する。
【００５９】
これにより、先頭フェーズｉｂｓより前方にＨ個の要素が追加され、収集開始時刻または所定時刻（収集開始時刻からＮフェーズ目）を基準にして、異なる時刻での組織の濃度値を含む複数の要素からなる第１の集合が１次元（時間軸方向）に配列され、（Ｈ＋Ｌ）×１の列ベクトルが構成される。
【００６０】
このシステム行列作成部２ｄは、組織に対応する時間濃度曲線を、第１の時刻（収集開始時刻または収集開始時刻からＮフェーズ目）を基準にして離散化し、異なる時刻での組織における造影剤の濃度に対応する値を含む複数の要素を１次元に配列した第１の集合を構成する第１の構成手段として機能する。
【００６１】
また、システム行列作成部２ｄは、動脈の時間濃度曲線Ａｉを時間軸方向に離散化することにより得た、異なる時刻での動脈の濃度値からなる集合の、時間軸方向の後方に、Ｈ個の要素を追加する。これにより、異なる時刻での動脈の濃度値を含む複数の要素からなる第２の集合が１次元（時間軸方向）に配列され、（Ｌ＋Ｈ）×１の列ベクトルが構成される。そしてシステム行列作成部２ｄは、１次元に配列された動脈の濃度値の列ベクトルをさらに２次元（空間方向）に配列した（Ｌ＋Ｈ）×（Ｌ＋Ｈ）の係数行列を構成する。
【００６２】
これにより、解析対象の先頭フェーズｉｂｓを基準にして、異なる時刻での動脈の濃度値を含む複数の要素からなる第２の集合が２次元（時間軸方向と空間方向）に配列され、（Ｌ＋Ｈ）×（Ｌ＋Ｈ）の係数行列が構成される。
【００６３】
このシステム行列作成部２ｄは、動脈に対応する時間濃度曲線を、第１の時刻（収集開始時刻または収集開始時刻からＮフェーズ目）より後方の第２の時刻（動脈の時間濃度曲線が上昇し始める時刻）を基準にして離散化し、異なる時刻での動脈における造影剤の濃度に対応する値を含む複数の要素を２次元に配列した第２の集合を構成する第２の構成手段としても機能する。
【００６４】
システム行列作成部２ｄは、動脈の時間濃度曲線とインパルス応答関数のデコンボリューションが組織の時間濃度曲線であることから、以上のように構成した、異なる時刻での組織の濃度値を含む複数の要素からなる列ベクトルＣ、異なる時刻での動脈の濃度値を含む複数の要素からなる係数行列Ａ、および未知数Ｘ（インパルス応答関数）を用いて、次式（１５）に示すようなＣ＝ＡＸのシステム行列を作成する。式（１５）において、Ｃは、（Ｈ＋Ｌ）×１ベクトルであり、Ａは、（Ｌ＋Ｈ）×（Ｌ＋Ｈ）行列であり、Ｘは、（Ｌ＋Ｈ）×１ベクトルである。ｉｂｓは、解析対象となる先頭フェーズを表わしている。
【数１５】
【００６５】
上記の式（１５）において、Ｈは、Ｌ×Ｌ係数行列に追加する要素の個数を表わしている。例えば、撮影（収集）開始時刻から先頭フェーズまでの範囲をフェーズ間隔で分割した個数、または、それより多くても少なくても構わない。追加要素の個数Ｈの数を少なくすれば、解析対象となる行列のサイズが小さく処理時間が短縮されるが、少なすぎると時間方向のずれに対する依存性が増大することとなる。従って、撮影開始時刻から先頭フェーズまでの範囲をフェーズ間隔で分割した個数の値を用いれば、時間方向のずれに対する依存性を小さくするのに十分な値となる。
【００６６】
組織の時間濃度曲線を時間軸方向に離散化した複数の要素のうち、追加分の要素値Ｐｉ（ｉ＝０，１，・・・Ｈ−１）には、０の要素を入力（追加）するが、部分的あるいは全部をＣｉとして、測定された組織の時間濃度曲線の値を用いるようにしても良い。この場合、実際にはＬ個以上のフェーズの画像データが撮影されている必要がある。
【００６７】
動脈の時間濃度曲線を時間軸方向に離散化した要素値ａi,jは、次式（１６）により表わされる。式（１６）において、ｉ＝０，１，・・・Ｌ−１であり、ｊ＝０，１，...Ｌ−１である。
【数１６】
【００６８】
上記の式（１６）において、追加分の要素値ｑｉ，ｊには、０の要素を入力（追加）するが、部分的あるいは全部を次式（１７）に示すようにして、測定された動脈の時間濃度曲線の値を用いるようにしても良い。この場合、実際にはＬ個以上のフェーズの画像データが撮影されている必要がある。
【数１７】
システム行列作成部２ｄは、作成した上記の式（１５）をデコンボリューション部２ｅに出力する。
【００６９】
デコンボリューション部２ｅは、システム行列作成部２ｄから出力された係数行列を入力し、SVD法やフーリエ変換法などの方法を用いてデコンボリューションを行う逆畳み込み積分手段として機能する。デコンボリューション部２ｅは、デコンボリューションにより得られた未知数Ｘから、上記の式（７）に示すようにしてインパルス応答関数Ｒｉを算出する。
【００７０】
デコンボリューション部２ｅは、上記の式（７）で算出されたインパルス応答関数Ｒｉから、ピークフェーズ（インパルス応答関数の最大値）の存在範囲を検索し、検索したピークフェーズからピーク範囲（幅）を算出する。
図３は、ピークフェーズを説明する図である。
【００７１】
図３に示すように、第１の実施形態において算出されるインパルス応答関数Ｒｉのピークは、従来のstandard SVD法に比べ、ほぼＨ個分後方のフェーズとなる。従って、デコンボリューション部２ｅは、ピークフェーズが後方にずれていることを考慮しなければならないため、局所血流動態に関する情報の算出の前に、予めピーク範囲を決定する必要がある。
【００７２】
ピーク範囲の算出方法として、例えば、先頭フェーズから時間軸方向にピークフェーズを検索し、検索したピークフェーズの前後で負の値の極小値をとるフェーズまでをピーク範囲とする。なお、負の値の極小値ではなく、０の値をとるフェーズまでをピーク範囲とするようにしてもよい。
【００７３】
デコンボリューション部２ｅは、検索したピークフェーズおよびピーク範囲に基づいて、局所血流動態に関する情報である、平均通過時間MTT（Mean Transit Time）、局所脳血流量CBF（Cerebral Blood Flow）、および局所脳血液量CBV（Cerebral Blood Volume）を、次式（１８）に示すようにして算出する。つまり、インパルス応答曲線のピークの高さから局所脳血流量CBFが求まり、インパルス応答曲線のピーク範囲の面積から局所脳血液量CBVが求まり、インパル応答曲線のピーク範囲から平均通過時間MTTが求まる。
［数１８］
CBF＝MAX(Ri）
CBV＝SUM(Ri)
MTT＝SUM(Ri)/MAX(Ri)・・・（１８）
【００７４】
このデコンボリューション部２ｅは、組織の伝達関数（インパルス応答関数）に基づいて、血流動態に関する情報を算出する血流情報算出手段としても機能する。
【００７５】
またデコンボリューション部２ｅは、算出した平均通過時間MTT、局所脳血流量CBF、および局所脳血液量CBVを次式（１９）に示すようにして補正する。
［数１９］
CBV=(100h/ρ)TRSUM(Ri)
CBF(r)=60(100h/ρ)MAX(Ri)
MTT2(r)=TR{SUM(Ri)/MAX(Ri)} ・・・（１９）
【００７６】
上記の式（１９）において、ρは、脳の比重を表わし、ｈは、ヘマトクリット値（一定量の血液中に含まれる赤血球の割合）を表わしている。なお、局所脳血流量CBFは、次式（２０）の関係を用いて算出（補正）するようにしてもよい。
［数２０］
CBF＝CBV/MTT2 ・・・（２０）
【００７７】
次に、図４のフローチャートを参照して、第１の実施形態における組織の血流動態に関する情報の算出処理について説明する。この処理では、特に、脳組織を例に挙げ、その血流動態に関する情報を算出する処理について説明する。
ステップＳ１において、時系列データ取得部２ａは、医用画像撮影装置が所定時間撮影して生成した時系列の撮影データを取得する。
【００７８】
ステップＳ２において、造影濃度曲線算出部２ｂは、ステップＳ１の処理で取得された時系列の撮影データから、各画素の画素値の変化を計測して造影濃度曲線を算出する。この造影濃度曲線を、組織の時間濃度曲線Ｃｉとする。
【００７９】
ステップＳ３において、動脈の時間濃度曲線算出部２ｃは、ステップＳ２の処理で算出された組織の時間濃度曲線Ｃｉから、CT画像やMRI画像の解析単位となる動脈領域を設定し、その動脈領域での時間濃度曲線ａiを算出する。このとき、動脈の時間濃度曲線算出部２ｃは、取得した撮影データがMRI画像の場合、ガンマ確率密度関数の定数倍の関数を動脈カーブモデルとして適用し、算出した動脈の時間濃度曲線ａｉのカーブフィッティングを行う。
ステップＳ４において、システム行列作成部２ｄは、解析対象となる先頭フェーズｉｂｓおよびフェーズ範囲ｉｂｓ≦ｉ≦Ｌ−１を決定する。
【００８０】
ステップＳ５において、システム行列作成部２ｄは、ステップＳ２の処理で算出された組織の時間濃度曲線Ｃｉのうち、ステップＳ４の処理で決定したフェーズ範囲ｉｂｓ≦ｉ≦Ｌ−１における組織の時間濃度曲線を時間軸方向に離散化し、離散化することにより得た異なる時刻での組織の濃度値からなる集合を１次元に配列して列ベクトルを構成する。
【００８１】
また、ステップＳ５において、システム行列作成部２ｄは、ステップＳ３の処理で算出された動脈の時間濃度曲線ａiのうち、ステップＳ４の処理で決定したフェーズ範囲ｉｂｓ≦ｉ≦Ｌ−１における動脈の時間濃度曲線を時間軸方向に離散化し、離散化することにより得た異なる時刻での動脈の濃度値からなる集合を２次元に配列して係数行列を構成する。
【００８２】
ステップＳ６において、システム行列作成部２ｄは、ステップＳ５の処理で構成された異なる時刻での組織の濃度値からなる列ベクトル（集合）の、時間軸方向の前方にＨ個の要素を追加する。これにより、解析対象の先頭フェーズｉｂｓが時間軸方向にＨ個分ずらされ、収集開始時刻または所定時刻を基準にして、異なる時刻での組織の濃度値を含む複数の要素からなる第１の集合が１次元に配列され、（Ｈ＋Ｌ）×１の列ベクトルＣが構成される。
【００８３】
また、ステップＳ６において、システム行列作成部２ｄは、ステップＳ５の処理で構成された異なる時刻での動脈の濃度値からなる係数行列（集合）の、時間軸方向の後方にＨ個の要素を追加する。これにより、動脈の造影濃度が上昇し始めるフェーズもしくは、それよりも前のフェーズを基準にして、異なる時刻での動脈の濃度値を含む複数の要素からなる第２の集合が２次元に配列され、（Ｌ＋Ｈ）×（Ｌ＋Ｈ）の係数行列Ａが構成される。
【００８４】
ステップＳ７において、システム行列作成部２ｄは、ステップＳ６の処理で得られた、異なる時刻での組織の濃度値を含む複数の要素からなる列ベクトルＣ、異なる時刻での動脈の濃度値を含む複数の要素からなる係数行列Ａ、および、未知数Ｘ（インパルス応答関数）を用いて、上記の式（１５）に示すシステム行列を作成する。
図５は、第１の実施形態に係る式（１５）のシステム行列を模式的に示す図である。
【００８５】
図５において、Ｃは、組織の時間濃度曲線Ｃｉを時間軸方向に離散化した（Ｈ＋Ｌ）×１の列ベクトルを模式的に表わしたものであり、Ａは、動脈の時間濃度曲線ａｉを時間軸方向に離散化した（Ｌ＋Ｈ）×（Ｌ＋Ｈ）の係数行列を模式的に表わしたものであり、Ｘは未知数（インパルス応答関数）を模式的に表わしたものである。
【００８６】
同図に示すように、組織の時間濃度曲線を構成する列ベクトルＣにおいて、前方にＨ個の要素ｐｉが追加され、動脈の時間濃度曲線を構成する係数行列Ａにおいて、後方にＨ個の要素ｑｉ，ｊが追加されている。追加要素の個数Ｈには、撮影開始時刻から先頭フェーズｉｂｓまでの範囲をフェーズ間隔で分割した個数の値を用いることが好ましい。追加分の要素値Ｐｉ（ｉ＝０，１，・・・Ｈ−１）と要素値ｑｉ，ｊには、０の要素が入力（追加）されるが、部分的あるいは全部に、実際の測定値を用いるようにしても良い。
【００８７】
このようなシステム行列を構成することで、従来のblock circulant SVD法のように係数行列を２倍にする必要がなくなり、デコンボリューションの高速処理化が可能になる。
【００８８】
図４の説明に戻る。ステップＳ８において、デコンボリューション部２ｅは、ステップＳ７の処理で作成されたシステム行列を、SVD法あるいはフーリエ変換法などの方法を用いてデコンボリューションを行い、ベクトルＸを求め、上記の式（７）に示すようにしてインパルス応答関数Ｒｉを算出する。
【００８９】
ステップＳ９において、デコンボリューション部２ｅは、ステップＳ８の処理で算出されたインパルス応答関数Ｒｉから、ピークフェーズを検索し、検索したピークフェーズからピーク範囲を算出する。
【００９０】
ステップＳ１０において、デコンボリューション部２ｅは、ステップＳ９の処理で算出されたピークフェーズおよびピーク範囲に基づいて、局所血流動態に関する情報である、平均通過時間MTT、局所脳血流量CBF、および局所脳血液量CBVを上記の式（１８）に従って算出する。
【００９１】
ステップＳ１１において、デコンボリューション部２ｅは、ステップＳ１０の処理で算出された平均通過時間MTT、局所脳血流量CBF、および局所脳血液量CBVを、上記の式（１９）に示すようにして補正する。
【００９２】
以上のように、血流動態解析装置１は、解析対象となる組織の時間濃度曲線を構成する列ベクトルの前方に所定個数の要素を追加し、動脈の時間濃度曲線を構成する係数行列の後方に所定個数の要素を追加し、組織の時間濃度曲線の先頭フェーズと、動脈の時間濃度曲線の先頭フェーズの基準位置をずらすことによって、処理時間を増大させることなく、かつ、時間方向のずれに対する依存性を小さくして、血管の血流動態に関する情報である、組織血流量、組織血液量、および平均通過時間を安定に算出することが可能となる。
【００９３】
これにより、病院側は、多数の患者の検査を実施することができるようになり、患者の検査待ち時間を短縮することができる。また、安定に解析結果を得ることができるため、再検査（再解析）といった事態を減らすことができ、患者への負担を軽減することができる。
【００９４】
また、血流動態解析装置１は、算出した平均通過時間MTT、局所脳血液量CBV、および局所脳血流量CBFに基づいて、対応する画像を表示部１ｈに表示させることができ、操作者は、血流動態の状態を容易に把握することが可能となる。
【００９５】
以上においては、対象組織として、脳組織を例に挙げ説明したが、第１の実施形態では、脳組織に限られるものではなく、肝臓、心臓、あるいは肺などの他の組織を対象とすることも勿論可能である。
【００９６】
また以上において、脳組織についてヘマトクリット値の補正および脳の比重による補正を行うようにした。そこで、例えば、肝臓や肺などの呼吸による動きがある組織、あるいは、心臓を対象組織とした場合については、組織の動き補正を行う必要がある。また動き補正だけでなく、ヘマトクリット値の補正、その他の補正を適宜行う必要がある。
【００９７】
（第２の実施形態）
図６は、本発明の第２の実施形態に係る流体解析装置の構成例を示す機能ブロック図である。
【００９８】
流体解析装置３は、例えば、MRI装置４に内蔵される。但し、ネットワークを介して流体解析装置３をMRI装置４或いはMRI装置４において収集された画像データを保管する画像処理装置や画像サーバに接続してもよい。
【００９９】
MRI装置４は、イメージング部４ａを備える。イメージング部４ａは、被検体の流体の非造影イメージング又は造影イメージングを実行することによって、複数フレーム分の非造影流体MR画像データ又は造影流体MR画像データを収集する機能を有する。流体としては、血管内の血流や脳脊髄液(CSF: cerebrospinal fluid)が挙げられる。MRI装置を用いた血管の撮影は、磁気共鳴血管撮影法(MRA: magnetic resonance angiography)と呼ばれる。
【０１００】
具体的には、イメージング部４ａは、静磁場磁石、傾斜磁場コイル、傾斜磁場電源、送信用高周波(RF: radio frequency)コイル、受信用RFコイル、送信器、受信器、シーケンサ及びコンピュータ等の構成要素を有する。そして、イメージング部４ａでは、コンピュータにおいて設定されたパルスシーケンス等の撮像条件に従って、静磁場磁石内の静磁場下にセットされた被検体に傾斜磁場コイルからは傾斜磁場が、送信用コイルからはRFパルスが、それぞれ印加される。そして、被検体から発生した核磁気共鳴(NMR: nuclear magnetic resonance)信号が受信用コイルで受信され、受信器においてデジタル化されたNMR信号に対する画像再構成処理がコンピュータにおいて実行されることによってMR画像データが生成される。また、傾斜磁場コイル及び送信用RFコイルには、それぞれコンピュータからシーケンサを通じて出力されるパルスシーケンスに従って傾斜磁場電源及び送信器から傾斜磁場パルス及びRFパルスが印加される。これにより、傾斜磁場コイル及び送信用RFコイルからは、所望の波形の傾斜磁場パルス及びRFパルスが被検体に印加される。
図７は、図６に示すMRI装置４のイメージング部４ａにおいて実行される非造影流体イメージング用の撮像条件の第１の例を示す図である。
【０１０１】
図７において横軸は時間を示す。図７は、非造影で周期性のないCSFを描出するためのパルスシーケンスの一例を示している。図７に示すように、ラベリングパルス（タグ付けパルスとも言う）をCSFまたは着目領域に印加し、続いて連続的に複数回に亘ってイメージングデータの収集(DATA ACQUISITION 1, DATA ACQUISITION 2, DATA ACQUISITION 3,...)を繰り返すダイナミックスキャンによって、異なる時刻t1, t2, t3, ...に対応する複数フレーム分の時系列の非造影流体画像データを収集することができる。
【０１０２】
ラベリングパルスとしては、t-SLIP (Time-SLIP: Time Spatial Labeling Inversion Pulse)、飽和(SAT: saturation)パルス（Rest grid pulseとも呼ばれる）、SPAMM (spatial modulation of magnetization)パルス及びダンテ(DANTE)パルスが知られている。同一種類のラベリングパルスを複数回印加したり、異種の複数のラベリングパルスを組み合わせて印加してもよい。
【０１０３】
例えば、複数の領域選択的90°SATパルスを組み合わせると、選択スラブ領域に放射状やストライプ状のパターンで飽和された領域を形成することができる。また、SPAMMパルスやダンテパルスの印加によってもストライプパターン、グリッドパターン（格子状のパターン）、放射状のパターン等の所望のパターンで飽和された領域を形成することができる。このため、パターンがCSFの流れに伴って移動する複数フレーム分のダイナミックCSF画像データを収集することができる。
【０１０４】
また、t-SLIPを用いたイメージングには、flow-in法とflow-out法がある。flow-in法は、t-SLIPを構成する領域選択反転回復(IR: inversion recovery)パルスを着目領域に印加し、着目領域の外部から着目領域に流入するラベリングされていないCSFを描出する方法である。一方、flow-out法は、t-SLIPを構成する領域非選択IRパルスを印加するとともにラベリング領域に領域選択IRパルスを印加し、ラベリング領域から着目領域に流入するラベリングされたCSFを描出する方法である。
【０１０５】
ラベリングパルスは、ECG (electro cardiogram)信号、呼吸センサによる同期信号又は脈波同期(PPG: peripheral pulse gating)信号等の生体信号に同期させて印加してもよい。
【０１０６】
一方、イメージングデータの収集用のシーケンスとしては、SSFP (steady state free precession)シーケンスやFASE (FastASE: fast asymmetric spin echoまたはfast advanced spin echo)シーケンス等の任意のシーケンスが使用できる。
図８は、図６に示すMRI装置４のイメージング部４ａにおいて実行される非造影流体イメージング用の撮像条件の第２の例を示す図である。
【０１０７】
図８において横軸は時間を示し、縦方向も時間を表している。図８は、IRパルスを用いて非造影で被検体の血流やCSF等の流体をイメージングするためのパルスシーケンスの一例を示している。すなわち、IRパルスの印加時刻からイメージングデータの収集タイミングまでの反転時間(TI: inversion time)をTI1, TI2, TI3, ...と徐々に変えながら、IRパルスの印加及びイメージングデータの収集(DATA ACQUISITION 1, DATA ACQUISITION 2, DATA ACQUISITION 3,...)が繰り返し実行される。
【０１０８】
IRパルスは単一のパルスに限らず、同時とみなせる十分に短い間隔で異種又は同種の複数のIRパルスを印加してもよい。領域選択的IRパルスは、上述のflow-in法の場合には撮像領域に印加され、flow-out法の場合には撮像領域の上流領域に印加される。そして、領域選択的IRパルス、周波数選択的IRパルス、非領域選択的IRパルス及び非周波数選択的IRパルスの組み合わせによって所望の対象物の磁化を倒すことによって撮像領域又は流体がラベリングされる。そして、TIに応じた距離だけ流体が移動したタイミングで撮像領域から収集したイメージングデータの画像再構成処理によって、撮像領域に流入する流体が選択的に描出された流体画像データを生成することができる。
【０１０９】
これにより、異なる複数のTIに対応する複数フレーム分の非造影流体画像データを収集することができる。すなわち、複数のTIに対応する複数フレーム分の画像データをTIの短い順に並べると、ダイナミック画像データのようにTIの増加に伴って流体が徐々に流れていく様子が描出される連続的な画像データを生成することができる。
【０１１０】
尚、IRパルスの印加及びイメージングデータの収集は、データ収集タイミングにおける流体の速度をより一定に近づけるために、ECG信号等の生体信号に同期させてもよい。生体信号を用いた同期イメージングは、血流のように周期性のある流体のイメージングに好適である。
【０１１１】
また、イメージングデータを収集するためのパルスシーケンスとしては、SE (spin echo)シーケンス、FSE (fast spin echo)シーケンス、FASEシーケンス、EPI (echo planar imaging) シーケンス等のパルスシーケンスが用いられる。
図９は、図６に示すMRI装置４のイメージング部４ａにおいて実行される非造影流体イメージング用の撮像条件の第３の例を示す図である。
【０１１２】
領域選択的IRパルスの印加及び図９に示すような血流を含む実線で示すイメージング領域からのイメージングデータの収集を繰り返し、かつTIを一定にしつつ領域選択的IRパルスの印加領域を点線で示すようにTAGGED REGION 1,TAGGED REGION 2, TAGGED REGION 3, ...に徐々に変えるイメージングによってもダイナミック画像データのように流体が徐々に流れていく様子が描出される連続的な画像データを生成することができる。すなわち、領域選択的IRパルスの印加領域に存在していた血流はラベリングされ、領域選択的IRパルスの印加領域から下流に向かってTIに応じた距離だけ流出したタイミングでイメージングデータが収集される。従って、ラベリングされた血流が選択的に描出された画像データを生成することができる。
【０１１３】
このため、領域選択的IRパルスの印加領域を変えて繰り返しイメージングデータを収集することによって、異なる複数のラベリング領域に対応し、ラベリングされた血流がイメージング領域の異なる位置に存在する複数フレーム分の非造影流体画像データを収集することができる。
【０１１４】
さらに、血管のより上流において血流が描出された画像データを各ラベリング領域に対応する画像データと最大値投影(MIP: maximum intensity projection)処理等の合成処理によって合成することにより、血流が徐々に流れていく様子が描出される連続的な画像データを生成することができる。
尚、図８及び図９においてIRパルスの代わりにSATパルスを用いてもよい。
【０１１５】
ここまでは、非造影イメージングの撮像条件について例を示したが、第１の実施形態と同様に公知の造影イメージング法によって造影画像データを収集することもできる。この場合、複数の時相に対応する時系列の複数フレーム分の造影画像データが収集される。
【０１１６】
一方、流体解析装置３は、イメージング部４ａにおいて収集された非造影流体MR画像データ又は造影流体MR画像データに基づいて、流体の動態を表す情報を取得する機能を有する。流体解析装置３は、例えば、コンピュータの演算装置に流体解析プログラムを読み込ませることにより、コンピュータをイメージデータ取得部３ａ、組織画素変化曲線算出部３ｂ、流体画素変化曲線算出部３ｃ、システム行列作成部３ｄ及びデコンボリューション部３ｅとして機能させたものである。但し、回路を用いて流体解析装置３の一部又は全部を構成してもよい。
【０１１７】
イメージデータ取得部３ａは、イメージング部４ａにおいて収集された複数フレーム分の連続的な流体画像データを取得して組織画素変化曲線算出部３ｂに与える機能を有する。
【０１１８】
組織画素変化曲線算出部３ｂは、複数フレーム分の連続的な流体画像データにおいて、組織に対応する複数の画素におけるフレーム間における画素値の変化から、組織内を通過するラベリングされた血流等の流体の濃度に応じた量の変化を求め、それぞれ組織の画素濃度変化曲線として算出する機能と、算出した複数の組織の画素濃度変化曲線を流体画素変化曲線算出部３ｃ及びシステム行列作成部３ｄに出力する機能を有する。
【０１１９】
尚、非造影イメージング法では、実際には造影剤が注入されない。このため、非造影イメージング法の場合には、血流等の流体の成分の量に関する濃度ではなく、ラベリングされた流体成分の流体全量に対する比に対応する論理的な濃度が求められる。
【０１２０】
例えば、時刻(t1, t2, t3, ...)、TI (TI1, TI2, TI3, ...)或いはラベリング領域(TAGGED REGION 1,TAGGED REGION 2, TAGGED REGION 3, ...)等のフレーム間で変えられるパラメータのi番目の値に対応する画素値をSiとし、初期値(t1, TI1, TAGGED REGION 1)等の基準となるパラメータの値に対応する画素値をS0とする。そうすると、第１の実施形態で示した式(１０)によって、組織における画素濃度に比例する量の変化を画素濃度変化曲線として算出することができる。或いは、フレーム間で変えられるパラメータに対応する組織の画素値S自体の変化を表す時間強度曲線(TIC: time intensity curve)を画素濃度変化曲線としてもよい。
画素濃度に比例する量及び画素値は、流体の流れに対応する量となる。従って、組織における画素濃度変化曲線は、組織における流体の動態情報を含んでいる。
【０１２１】
流体画素変化曲線算出部３ｃは、組織の複数の画素における画素濃度変化曲線に基づいて、動脈等の流体領域内を循環する流体の画素濃度変化曲線を算出する機能と、算出した流体の画素濃度変化曲線をシステム行列作成部３ｄに出力する機能を有する。流体の画素濃度変化曲線は、流体領域内を循環する流体に対応する画素値の濃度に応じた量のフレーム間の変化を表す曲線である。
【０１２２】
流体の画素濃度変化曲線は、第１の実施形態と同様な方法で算出することができる。このため、第１の実施形態と同様に組織の画素濃度曲線は、複数の画素についてそれぞれ求められるのに対して、流体の画素濃度曲線は流体の画素値の濃度に応じた量を代表する値を用いて１つの曲線として求められる。従って、流体部位として複数の画素を含む関心領域(ROI: region of interest)が設定された場合には、第１の実施形態と同様に画素値の濃度に応じた量を代表する値が算出される。代表値としては、画素間の平均値を用いる方法や式（１１）に示すようなフィッティングパラメータを有するガンマ確率密度関数の定数倍の関数を用いたカーブフィッティングにより算出された値を用いる方法がある。流体領域に対応する画素についても、第１の実施形態と同様に自動的あるいはユーザからの指定によって特定することができる。
【０１２３】
システム行列作成部３ｄは、組織の各画素に対応する画素濃度変化曲線に基づいて、時刻(t1, t2, t3, ...)、TI (TI1, TI2, TI3, ...)或いはラベリング領域(TAGGED REGION 1,TAGGED REGION 2, TAGGED REGION 3, ...)等のパラメータの値ごとの組織の画素値の濃度に応じた量を要素とする１次元の列ベクトルＣを作成する機能と、流体の画素濃度変化曲線に基づいて、パラメータの値ごとの流体の画素値の濃度に応じた量を要素とする１次元の列ベクトルを作成する機能とを有する。加えて、システム行列作成部３ｄは、組織に対応する列ベクトルＣ、流体に対応する列ベクトルを２次元に配列した係数行列Ａ及び未知数Ｘを用いて、Ｃ＝ＡＸのシステム行列を作成する機能と、作成したシステム行列をデコンボリューション部３ｅに出力する機能とを有する。
【０１２４】
組織に対応する列ベクトルＣ、流体に対応する列ベクトル、係数行列Ａ及びシステム行列Ｃ＝ＡＸは、それぞれ第１の実施形態における方法と同様な方法で作成することができる。換言すれば以下のような説明になる。
図１０は、図６に示すシステム行列作成部３ｄにおける解析対象となる要素の作成方法を説明する図である。
【０１２５】
図１０において横軸は、時刻(t1, t2, t3, ...)、TI (TI1, TI2, TI3, ...)或いはラベリング領域(TAGGED REGION 1,TAGGED REGION 2, TAGGED REGION 3, ...)等のパラメータ値を示し、縦軸は画素値の濃度を示す。
【０１２６】
図１０に示すように、組織の各画素に対応する画素濃度曲線C1i, C2i, C3i, ...及び流体の画素濃度曲線aiが得られている場合に、複数の組織の画素濃度変化曲線C1i, C2i, C3i, ...がそれぞれ少なくとも立ち上がる前となる第１の点及び流体の画素濃度変化曲線aiが少なくとも立ち上がる前となる第２の点が決定される。例えば、第１の点は、時刻(t1, t2, t3, ...)、TI (TI1, TI2, TI3, ...)或いはラベリング領域(TAGGED REGION 1,TAGGED REGION 2, TAGGED REGION 3, ...)等のイメージングにおいて変えられたパラメータの初期値(t1, TI1, TAGGED REGION 1)に決定される。また、第２の点は、第１の点よりも順番が後の値に決定される。
【０１２７】
さらに、組織の画素濃度曲線C1i, C2i, C3i, ...に対応する解析対象となる要素と流体の画素濃度曲線aiに対応する解析対象となる要素が同じ数になるように、第１の点及び第２の点を開始点とする同じ広さの解析対象範囲RANGE1, RANGE2が組織の画素濃度曲線C1i, C2i, C3i, ...及び流体の画素濃度曲線aiに対してそれぞれ設定される。尚、組織の画素濃度曲線C1i, C2i, C3i, ...の解析対象範囲RANGE1のみに含まれるパラメータ値における前方の一部又は全部の要素をゼロとしてもよい。一方、流体の画素濃度曲線aiの解析対象範囲RANGE2のみに含まれるパラメータ値における後方の一部又は全部の要素をゼロとしてもよい。
【０１２８】
この結果、第１の点における組織の画素値の濃度に応じた量又はゼロが組織の画素濃度曲線C1i, C2i, C3i, ...に対応する解析対象となる最初の要素となり、第２の点における流体の画素値の濃度に応じた量が流体の画素濃度曲線aiに対応する解析対象となる最初の要素となる。
【０１２９】
このように作成された同じ数の組織の画素濃度曲線に対応する要素及び流体の画素濃度曲線に対応する要素を用いて式（１５）に示すようなシステム行列を作成することができる。
【０１３０】
デコンボリューション部３ｅは、第１の実施形態におけるデコンボリューション部２ｅが有する機能と同様の機能を有する。すなわち、デコンボリューション部３ｅは、組織に対応する列ベクトルＣと流体に対応する列ベクトルを２次元に配列した係数行列Ａの逆畳み込み積分によって未知数Ｘを算出する機能と、未知数Ｘから得られるインパルス応答関数に基づいて流体の平均通過時間MTT、局所脳血流量CBF、局所脳血液量CBV、局所CSF流量、局所CSF量等の流体動態に関する情報を算出する機能を有する。
【０１３１】
また、必要に応じてデコンボリューション部３ｅは、第１の実施形態におけるデコンボリューション部２ｅと同様に、流体動態に関する情報を式（１９）に示すように補正するように構成される。
次に、流体解析装置３の動作及び作用について説明する。
【０１３２】
図１１は、図６に示す流体解析装置３におけるデータ処理の流れを示すフローチャートである。尚、第１の実施形態において図４に示されている詳細な処理と同様な処理については省略する。
まず予め、MRI装置４のイメージング部４ａは非造影イメージング又は造影イメージングを実行し、複数フレーム分のMR流体画像データを収集する。
【０１３３】
そして、ステップＳ２１において、流体解析装置３のイメージデータ取得部３ａは、イメージング部４ａから複数フレーム分の流体画像データを取得して組織画素変化曲線算出部３ｂに与える。
【０１３４】
次に、ステップＳ２２において、組織画素変化曲線算出部３ｂは、複数フレーム分の流体画像データの組織に対応する複数の画素におけるフレーム間における画素値の濃度に応じた量の変化を求めることによって組織の画素濃度変化曲線を算出する。
【０１３５】
次に、ステップＳ２３において、流体画素変化曲線算出部３ｃは、組織の複数の画素における画素濃度変化曲線に基づいて、流体の画素濃度変化曲線を算出する。
【０１３６】
次に、ステップＳ２４において、システム行列作成部３ｄは、組織の画素濃度変化曲線及び流体の画素濃度変化曲線に基づいて、組織の画素濃度に対応する列ベクトルをＣ、流体の画素濃度に対応する列ベクトルを２次元に配列した係数行列をＡ、未知数をＸとするＣ＝ＡＸのシステム行列を作成する。このとき、組織に対応する列ベクトルＣの最初の要素は、組織の画素濃度変化曲線の横軸上の第１のパラメータ値に対応する値とされる。一方、流体に対応する列ベクトルの最初の要素は、流体の画素濃度変化曲線の横軸上の第１のパラメータ値よりも後方の第２のパラメータ値に対応する値とされる。また、組織に対応する列ベクトルＣの要素数と流体に対応する列ベクトルの要素数が同じになるように各列ベクトルが作成される。そのために、列ベクトルを構成する要素として画素濃度の代わりにゼロを用いてもよい。
【０１３７】
次に、ステップＳ２５において、デコンボリューション部３ｅは、システム行列の未知数Ｘを算出し、未知数Ｘから得られるインパルス応答関数に基づいて流体動態に関する情報を算出する。この流体動態に関する情報は、式（１９）によって補正してもよい。
【０１３８】
このように、流体解析装置３では、複数の時相に対応する時系列の造影画像データに限らず、撮像パラメータを変えて収集された連続的な複数フレーム分の非造影画像データに基づいて、造影画像データに対する処理と同様な処理で血流やCSF等の流体の動態に関する情報を求めることができる。
【０１３９】
尚、第２の実施形態では、流体解析装置３がMRI装置で収集された画像データに基づいて流体動態に関する情報を取得する例について説明したが、Ｘ線CT装置、超音波診断装置又はPET装置等の医用画像撮影装置で収集された造影画像データ又は非造影画像データに基づいて流体解析装置３が流体動態に関する情報を取得するように構成することもできる。
【０１４０】
（他の実施形態）
以上、特定の実施形態について記載したが、記載された実施形態は一例に過ぎず、発明の範囲を限定するものではない。ここに記載された新規な方法及び装置は、様々な他の様式で具現化することができる。また、ここに記載された方法及び装置の様式において、発明の要旨から逸脱しない範囲で、種々の省略、置換及び変更を行うことができる。添付された請求の範囲及びその均等物は、発明の範囲及び要旨に包含されているものとして、そのような種々の様式及び変形例を含んでいる。
【符号の説明】
【０１４１】
１ 血流動態解析装置
２ 局所血流解析アプリケーション
２ａ 局所血流解析アプリケーション
２ｂ 造影濃度曲線算出部
２ｃ 動脈の時間濃度曲線算出部
２ｄ システム行列作成部
２ｅ デコンボリューション部
３ 流体解析装置
３ａ イメージデータ取得部
３ｂ 組織画素変化曲線算出部
３ｃ 流体画素変化曲線算出部
３ｄ システム行列作成部
３ｅ デコンボリューション部
４ MRI装置
４ａ イメージング部
【技術分野】
【０００１】
本発明の実施形態は、血流動態解析装置、血流動態解析プログラム、流体解析装置及び流体解析プログラムに関する。
【背景技術】
【０００２】
従来、Ｘ線CT（Computed Tomography）装置やMRI（Magnetic Resonance Imaging）装置では、静脈から造影剤を投与し、時系列画像データを収集して、その画像を元に解析し、組織内の血流動態に関する情報を得る検査が行われている。これはパーヒュージョン検査と呼ばれ、造影剤の撮影断面への集中の度合いが画像の濃度変化として収集できることを利用したものである。
【０００３】
例えば、肺循環や造影剤投与のばらつきをなくすために、組織に流入する動脈の時間濃度曲線（Time Concentration Curve：TCC）：ａｉを入力関数として、測定した組織の時間濃度曲線Ｃｉとのデコンボリューション（Deconvolution：逆畳み込み積分）を行い、組織固有のインパルス応答関数（伝達関数）Ｒｉ（residue function）を求め、得られたインパルス応答関数から、血流動態を定量的に表す指標である血流量（脳の場合、CBF：Cerebral Blood Flow）、平均通過時間：MTT（Mean Transit Time）、および、血液量（脳の場合、CBV：Cerebral Blood Volume）などのパラメータを算出している。この解析方法は、デコンボリューション法と呼ばれている。
【０００４】
尚、時間濃度曲線は、造影剤濃度に比例する量の時間変化を表す曲線である。Ｘ線CT装置で収集された造影画像データから求められる造影剤濃度に比例する量の時間変化曲線は、TDC (Time Density Curve)とも呼ばれるが、ここではTCCと表記する。
【０００５】
例えば、非特許文献１では、デコンボリューションに、standard SVD（Singular Value Decomposition）法やblock circulant SVD法を適用している。
【０００６】
ここで、standard SVD法について説明する。まず、収集された画像データにおいて、画素ｊ（位置ｒｊ）の画素値の変化を計測して組織の時間濃度曲線Ｃｉ（ｒｊ）を算出する。また、算出した組織の時間濃度曲線Ｃｉから、動脈の時間濃度曲線Ａｉを次式（１）に示すようにして算出する。式（１）において、ｉは、フェーズ番号を表わし、ΩＡは、動脈に相当する画素の集合を表わしている。
【数１】
【０００７】
組織の時間濃度曲線Ｃｉが、動脈の時間濃度曲線Ａｉとインパルス応答関数Ｘのデコンボリューション（畳み込み積分）であるとすると、次式（２）に示すようなシステム方程式で表わすことができる。式（２）において、ＴＲは、画像撮影の時間間隔（繰り返し時間 TR: repetition time)を表わし、Ｒｊは、インパルス応答関数を表わしている。
【数２】
【０００８】
例えば、動脈の時間濃度曲線ａｉおよび組織の時間濃度曲線Ｃｉを時間軸方向に離散化したもの、つまり、ｉｂｓ≦ｉ≦Ｌ−１のフェーズ範囲をＬ個の等区間に分割し、Ｌ個のフェーズのデータを解析に用いるとすると、上記の式（２）は、行列形式でＣ＝ＡＸとなり、次式（３）に示すように表わすことができる。式（３）において、Ｃは、Ｌ×１の列ベクトルであり、Ａは、Ｌ×Ｌの係数行列であり、Ｘは、Ｌ×１の列ベクトルである。
【数３】
【０００９】
上記の式（３）において、ｉｂｓは、解析対象となるフェーズ範囲の先頭フェーズを表わしており、動脈および組織の造影濃度が上昇し始めるフェーズ、もしくは、それよりも前のフェーズを用いる。また、係数行列の各要素値ａi,jは、次式（４）により表わされる。式（４）において、ｉ＝０，１，・・・Ｌ−１であり、ｊ＝０，１，...Ｌ−１である。
【数４】
そして、上記の式（３）の近似解を、特異値分解により求める。すなわち、行列Ａを次式（５）に示すように３つの行列に分解することができる。
［数５］
Ａ＝ＵＤＶＴ ・・・（５）
【００１０】
上記の式（５）において、Ｕは、ｍ行ｋ列の正規直交行列であり、Ｄは、ｋ行ｋ列の対角行列であり、Ｖは、ｎ行ｋ列の正規直交行列である。上記の式（３）に用いる場合、ｍ＝ｎ＝ｋ＝Ｌである。対角行列とは、対角成分以外の成分がすべて０である行列であり、正規直交行列とは、対角成分が１である対角行列である。
【００１１】
ここで、Ｕ＝（ｕ１，ｕ２，...ｕｋ）、Ｖ＝（ｖ１，ｖ２，...ｖｋ）、Ｄ＝diag（ｗ０，ｗ１，...ｗｋ）であるとき、ｗｋを降順に並び替え、ｕｋ、ｖｋもその順序に対応するように入れ替えを行っておく。行列Ａの特異値分解の結果を用いて、Ｘの最小２乗解である最小ノルム推定値を次式（６）により算出することができる。
【数６】
【００１２】
上記の式（６）において、Ｗは、Ｄの対角要素（特異値）の逆数からなる対角行列であって、Ｗ＝diag（１／ｗｏ，１／ｗ１，...１／ｗＧ−１，０，０，...）である。ここで、最大特異値のＰｓｖｄ倍に満たない特異値に対しては、逆数ではなく０を対角要素とする。従って、Ｇ個の対角要素だけが０でない値を持つようになる。Ｐｓｖdは、正則化の強さを左右するパラメータであり、大きい値ほど強い正則化が行われる。
そして、ノルム推定値が求められたら、次式（７）に示すようにしてインパルス応答関数Ｒｉを算出することができる。
【数７】
上記の式（７）で算出されたインパルス応答関数Ｒｉから、組織血流量、組織血液量、および平均通過時間を算出することができる。
【００１３】
ところで、血流は、組織の毛細血管へ流れるが、造影剤が組織に到達する時間は、組織の場所によって異なり、設定した動脈領域から求めた動脈の時間濃度曲線と各画素の組織の時間濃度曲線は、厳密には時間軸が一致しない。従って、極端な場合、動脈の造影濃度が上昇する前に（動脈の時間濃度曲線が立ち上がる前に）、組織の造影濃度が上昇してしまう（組織の時間濃度曲線が先に立ち上がってしまう）こととなる。
【００１４】
上述したstandard SVD法では、このような時間方向のずれ（ディレイ）がないことを仮定としている。しかしながら、動脈の時間濃度曲線と組織の時間濃度曲線の間には、必ず時間方向のずれが存在するため、ディレイ時間の値によってデコンボリューション結果が変わり、算出される組織血流量の値が影響を受けてしまうこととなる。
【００１５】
そこで、block circulant SVD法では、時間方向のずれに対する依存性を軽減するために、次式（８）に示すようなブロック循環行列であるシステム行列を用いるようにする。
【数８】
【００１６】
このblock circulant SVD法において、基本的な処理は、standard SVD法と同様であり、行列形式Ｃ＝ＡＸであるが、Ａが、２Ｌ×２Ｌの係数行列となっている。そのため、動脈の時間濃度曲線を構成する係数行列Ａの各要素値の後方に０が追加され、インパルス応答関数を構成する列ベクトルＸの各要素値の後方に０が追加され、組織の時間濃度曲線を構成する列ベクトルＣの各要素値の後方に０が追加され、使用する測定データの２倍の要素数となるように拡張される。特に、係数行列Ａは、行および列がいずれも２倍、もしくはそれ以上に拡張され、さらに列ベクトルの要素が周期的になるように配列し直されている。
【００１７】
このように、システム行列を上記の式（８）に示すように構成することで、時間方向のずれが存在した場合、インパルス応答関数が時間軸方向にシフトし、その形状が相似形を保つようになる。従って、インパルス応答関数を時間軸方向にシフトさせることによって、時間方向のずれに対する依存性を小さくした、組織血流量、組織血液量、および平均通過時間を算出することができる。
【００１８】
しかしながら、システム行列において、動脈の時間濃度曲線を構成する係数行列Ａが２×２倍、組織の時間濃度曲線を構成する列ベクトルが２倍となることから、解析処理時間は８倍以上にも増大してしまい、実用的ではない。
【００１９】
そこで、非特許文献２では、解析処理時間を増大させずに、standard SVD法のシステム行列を応用している。すなわち、非特許文献２で提案されている方法では、システム行列として上記の式（４）を用い、既知ベクトルには、先頭フェーズからオフセット位置までの組織の時間濃度曲線の要素をノイズ（信号のひずみ）として計算上から除外している。ここでは、上記の式（３）において、ベクトルＣからオフセット位置までの先頭数フェーズが除かれることから、後方には、同数の０の要素が追加されることになる。
このように、非特許文献２の技術では、システム行列を拡張せずに、時間方向のずれに対する依存性を小さくするようにしている。
【先行技術文献】
【非特許文献】
【００２０】
【非特許文献１】Ona Wu et al,"Tracer Arrival Timing-Insensitive Technique for Estimating Flow in MR Perfusion-Weighted Imaging Using Singular Value Decomposition With a Block-Circulant Deconvolution Matrix",Magnetic Resonance in Medicine 50:164-174(2003)
【非特許文献２】M.R.Smith, H.Lu, S.Trochet, and R.Frayne "Removing the Effect of SVD Algorithmic Artifacts Present in Quantitative MR Perfusion Studies", Magnetic Resonance in Medicine 51:631-634(2004)
【発明の概要】
【発明が解決しようとする課題】
【００２１】
非特許文献１の技術において、standard SVD法では、時間方向のずれを考慮していないため、デコンボリューション結果の正確性に欠ける。一方、block circulant SVD法では、時間方向のずれを考慮しているものの、解析処理時間が膨大であり、いずれも実用的ではない課題があった。
【００２２】
特に、Ｘ線CT装置やMRI装置の運用では、通常、多数の患者の検査を実施しなければならない。また、画像撮影後、直ちに解析処理を行って、技師が、処理結果を見て撮影に不備があるか否かを確認しなければならない。従って、解析処理時間が膨大になると、実施可能な検査数が減少してしまうとともに、患者にとっては待ち時間が増え、病院にとっては患者獲得の減少になり、双方に不利益をもたらすこととなる。
【００２３】
そこで、非特許文献２の技術では、解析処理時間を増大させずに処理を行っているが、先頭フェーズからオフセット位置までの組織の時間濃度曲線の要素をノイズ（信号のひずみ）として計算上から除外してしまうことから、オフセット位置に結果が依存してしまうという新たな問題が生じてしまう。しかしながら、オフセット位置を正しく同定することは難しいため、結果的に血流動態に関する情報を安定に求めることができないという課題があった。
【００２４】
本発明はこのような状況に鑑みてなされたものであり、その目的は、処理時間を増大させることなく流体の動態に関する情報を安定に算出することが可能な血流動態解析装置、血流動態解析プログラム、流体解析装置及び流体解析プログラムを提供することである。
【課題を解決するための手段】
【００２５】
本発明の実施形態に係る血流動態解析装置は、組織の時間濃度曲線算出手段、動脈の時間濃度曲線算出手段、第１の構成手段、第２の構成手段、逆畳み込み積分手段及び血流情報算出手段を備える。組織の時間濃度曲線算出手段は、医用画像撮影装置により得られた複数の時相における被検体の造影画像データに基づいて、前記複数の時相間における各画素の画素値の変化を計測して組織における造影剤の濃度に対応する量の時間変化を表す時間濃度曲線を算出する。動脈の時間濃度曲線算出手段は、前記組織に対応する前記時間濃度曲線から動脈領域を設定し、動脈に対応する時間濃度曲線を算出する。第１の構成手段は、前記組織に対応する前記時間濃度曲線を、第１の時刻を基準にして離散化し、異なる時刻での前記組織における前記造影剤の濃度に対応する値を含む複数の要素を１次元に配列した第１の集合を構成する。第２の構成手段は、前記動脈に対応する前記時間濃度曲線を、前記第１の時刻より後方の第２の時刻を基準にして離散化し、異なる時刻での前記動脈における前記造影剤の濃度に対応する値を含む複数の要素を２次元に配列した第２の集合を構成する。逆畳み込み積分手段は、前記第１の集合と前記第２の集合の逆畳み込み積分を行い、前記組織の伝達関数を算出する。血流情報算出手段は、前記伝達関数に基づいて、血流動態に関する情報を算出する。
【００２６】
また、本発明の実施形態に係る血流動態解析プログラムは、コンピュータを、上述の組織の時間濃度曲線算出手段、動脈の時間濃度曲線算出手段、第１の構成手段、第２の構成手段、逆畳み込み積分手段及び血流情報算出手段として機能させる。
【００２７】
また、本発明の実施形態に係る流体解析装置は、濃度変化算出手段及び流体情報算出手段を備える。濃度変化算出手段は、被検体の複数のフレームの画像データに基づいて、組織に対応する画素において画素値の濃度に応じた量の前記複数のフレーム間における第１の変化を算出する一方、流体に対応する画素値の濃度に応じた量の前記複数のフレーム間における第２の変化を算出する。流体情報算出手段は、前記第１の変化に基づいて前記第１の変化の第１の点に対応する要素を最初の要素とする解析対象となる数の第１の複数の要素を作成する一方、前記第２の変化に基づいて前記第２の変化の前記第１の点よりも後の第２の点に対応する要素を最初の要素とする前記解析対象となる数の第２の複数の要素を作成し、前記第１の複数の要素及び前記第２の複数の要素に基づいて逆畳み込み積分処理を伴って流体の動態に関する情報を算出する。
また、本発明の実施形態に係る流体解析プログラムは、コンピュータを、上述の濃度変化算出手段及び流体情報算出手段として機能させる。
【図面の簡単な説明】
【００２８】
【図１】本発明の第１の実施形態に係る血流動態解析装置の構成例を示す図。
【図２】局所血流解析アプリケーションを実行するCPUの機能の一例を示すブロック図。
【図３】ピークフェーズを説明する図。
【図４】組織の血流動態に関する情報の算出処理を説明するフローチャート。
【図５】第１の実施形態に係るシステム行列を模式的に示す図。
【図６】本発明の第２の実施形態に係る流体解析装置の構成例を示す機能ブロック図。
【図７】図６に示すMRI装置のイメージング部において実行される非造影流体イメージング用の撮像条件の第１の例を示す図。
【図８】図６に示すMRI装置のイメージング部において実行される非造影流体イメージング用の撮像条件の第２の例を示す図。
【図９】図６に示すMRI装置のイメージング部において実行される非造影流体イメージング用の撮像条件の第３の例を示す図。
【図１０】図６に示すシステム行列作成部３ｄにおける解析対象となる要素の作成方法を説明する図。
【図１１】図６に示す流体解析装置におけるデータ処理の流れを示すフローチャート。
【発明を実施するための形態】
【００２９】
以下、本発明の実施の形態について図面を参照して詳細に説明する。
【００３０】
（第１の実施形態）
図１は、本発明の第１の実施形態に係る血流動態解析装置１の構成例を示す図である。
【００３１】
血流動態解析装置１は、CPU(Central Processing Unit)１ａ、ROM（Read Only Memory）１ｂ、RAM（Random Access Memory）１ｃ、および入出力インターフェイス１ｄが、バス１ｅを介して接続されている。入出力インターフェイス１ｄには、記憶部１ｆ、通信部１ｇ、表示部１ｈ、入力部１ｉ、およびリムーバブルメディア１ｊが接続されている。
【００３２】
CPU１ａは、入力部１ｉからの入力信号に基づいて血流動態解析装置１を起動するためのブートプログラムをROM１ｂから読み出して実行し、記憶部１ｆに格納されている各種オペレーティングシステムを読み出す。またCPU１ａは、入力部１ｉからの入力信号に基づいて各種の制御を行ったり、ROM１ｂや記憶部１ｆに記憶されたプログラムおよびデータを読み出してRAM１ｃにロードしたり、あるいはRAM１ｃから読み出されたプログラムのコマンドに基づいて、データ演算または加工などの一連の処理を実行する。
【００３３】
記憶部１ｆは、半導体メモリや磁気ディスクなどで構成されており、CPU１ａで実行されるプログラムやデータを記憶する。記憶部１ｆには、CPU１ａが実行するプログラムとして、例えば、局所血流動態に関する情報を解析するためのアプリケーションが用意される。以下、局所血流動態に関する情報を解析するこのアプリケーションを局所血流解析アプリケーションという。
【００３４】
通信部１ｇは、LAN（Local Area Network）カードやモデムなどで構成されており、血流動態解析装置１をLANやインターネットといった通信媒体に接続することを可能にする。すなわち通信部１ｇは、通信媒体から受信したデータを、入出力インターフェイス１ｄおよびバス１ｅを介してCPU１ａに送信し、CPU１ａからバス１ｅおよび入出力インターフェイス１ｄを介して受信したデータを、通信媒体に送信する。
【００３５】
表示部１ｈは、例えば液晶ディスプレイである。CPU１ａからバス１ｅおよび入出力インターフェイス１ｄを介して受信した信号に基づいて、CPU１ａの処理結果などを表示する。
【００３６】
入力部１ｉは、血流動態解析装置１の操作者が各種の操作を入力するキーボードやマウスなどの入力デバイスにより構成されている。操作者の操作に基づいて入力信号を生成し、入出力インターフェイス１ｄおよびバス１ｅを介してCPU１ａに送信する。
【００３７】
リムーバブルメディア１ｊは、例えば光ディスクやフレキシブルディスクなどである。図示せぬディスクドライブによって読み出されたデータが、入出力インターフェイス１ｄおよびバス１ｅを介してCPU１ａに送信され、CPU１ａからバス１ｅおよび入出力インターフェイス１ｄを介して受信したデータが、ディスクドライブによって書き込まれる。
【００３８】
第１の実施形態における血流動態解析装置１は、Ｘ線CT装置やMRI装置などの医用画像撮影装置（いわゆる医用モダリティ）により、ダイナミックス撮影により収集された画像データから血流動態を解析する装置である。このため、血流動態解析装置１は、かかる画像データを入手できる環境にあればよく、医用画像撮影装置と一体に構成されていてもよいし、医用画像撮影装置とは別体で構成されていてもよい。別体で構成される場合には、収集された画像データはリムーバブルメディア１ｊで、または通信部１ｇを介して図示せぬ医用画像撮影装置から血流動態解析装置１に送られる。
【００３９】
また、血流動態解析装置１は、血流動態解析結果を、通信部１ｇを介して、PACS(Picture Archiving and Communication System)に転送し、そこに保存させることができる。これによって、特定の画像データの血流動態解析結果の検索や閲覧などを容易に行うことが可能となる。
【００４０】
医用画像撮影装置には、Ｘ線CT装置やMRI装置のほか、超音波診断装置やPET（Positron Emission Tomography）装置などが含まれる。
図２は、局所血流解析アプリケーションを実行するCPU１ａの機能の一例を示すブロック図である。
【００４１】
図２に示すように、局所血流解析アプリケーションを実行するCPU１ａは、時系列データ取得部２ａ、造影濃度曲線算出部２ｂ、動脈の時間濃度曲線算出部２ｃ、システム行列作成部２ｄ、およびデコンボリューション部２ｅとして機能する。
【００４２】
時系列データ取得部２ａは、造影剤を静脈注入（ボーラス注入）した被検者について、医用画像撮影装置が所定時間撮影して生成した時系列の撮影データ（多時相の画像データ）を取得し、取得した時系列の撮影データを造影濃度曲線算出部２ｂに出力する。
造影濃度曲線算出部２ｂは、時系列データ取得部２ａから出力される時系列の撮影データを入力し、各画素の画素値の変化を計測して造影濃度曲線を算出する。
【００４３】
つまり、被検者に造影剤を比較的短時間で静脈注入すると、造影剤は静脈を通り心臓に戻り、肺を循環して再び心臓に戻り、動脈へ流入し、各臓器の組織に到達し、再び静脈を通り心臓へ戻る。ここで、例えば、医用画像撮影装置が脳の領域のCT画像やMRI画像を時系列で撮影していた場合、造影剤濃度の変化に応じた各画素の画素値の変化を計測することができる。造影濃度曲線算出部２ｂは、この計測した各画素の画素値の変化に基づいて、造影濃度曲線を算出することができる。
【００４４】
例えば、CT画像の場合、次式（９）に示すように、各時相の画素値Ｓｉから造影剤到着前の画素値Ｓ０を引き算することにより、造影剤濃度に比例する量（TDC：Time Density Curve）が算出される。式（９）において、ｉは、撮影時相を表している。
［数９］
Ｃｉ＝Ｓｉ−Ｓ０ ・・・（９）
【００４５】
また例えば、MRI画像の場合、次式（１０）に示すように、各時相の画素値Ｓｉと造影剤到着前の画素値Ｓ０の比の対数（ΔR2＊）から、造影剤濃度に比例する量（TCC：Time Concentration Curve）が算出される。式（１０）において、ｉは、撮影時相を表し、TEは、エコー時間を表している。
［数１０］
ΔＲ２＊＝−log（Ｓｉ／Ｓ０）／TE ・・・（１０）
【００４６】
造影濃度曲線算出部２ｂは、上記の式（９）または式（１０）で算出した造影濃度曲線を、組織の時間濃度曲線Ｃｉとし、算出値を動脈の時間濃度曲線算出部２ｃ、および、システム行列作成部２ｄに出力する。
【００４７】
この造影濃度曲線算出部２ｂは、複数の時相における被検体の造影画像データに基づいて、複数の時相間における各画素の画素値の変化を計測して組織における造影剤の濃度に対応する量の時間変化を表す時間濃度曲線を算出する組織の時間濃度曲線算出手段として機能する。
【００４８】
動脈の時間濃度曲線算出部２ｃは、造影濃度曲線算出部２ｂから出力される組織の時間濃度曲線Ｃｉを入力し、その曲線から、CT画像やMRI画像の解析単位となるボクセルでの動脈の時間濃度曲線を算出する。つまり、CT画像やMRI画像の解析単位となるボクセルの大きさは、０．５ｍｍ乃至３ｍｍ程度の大きさであり、組織の中の主要な動脈であれば、血管にほぼ包含されるボクセルが存在する。
【００４９】
動脈の時間濃度曲線算出部２ｃは、このようなボクセルを自動的あるいはユーザからの指定によって特定することで、動脈領域を設定し、その動脈領域での時間濃度曲線（ここでは、ｉ番目の撮影時刻での濃度値ａi）を算出する動脈の時間濃度曲線算出手段として機能する。なお、複数の画素が動脈領域として特定された場合には、例えば、その平均値をａiとする。
【００５０】
或いは、動脈の時間濃度曲線算出部２ｃは、必要に応じて、次式（１１）によるガンマ確率密度関数の定数倍の関数（ガンマ確率密度関数に比例する関数）を動脈の時間濃度曲線のカーブモデル（以下、「動脈カーブモデル」と称する）として適用し、算出した動脈の時間濃度曲線ａｉのカーブフィッティング（動脈カーブモデルへの当てはめ）を行う。通常、MRI画像ではカーブフィッティングを行う場合が多いが、CT画像ではカーブフィッティングを行わない場合が多い。
［数１１］
φ(t; t0, α, β, K)=K(t-t0)αexp{-β-1(t-t0)} ・・・（１１）
【００５１】
ここで言うガンマ確率密度関数とは、時間βに１回の割合で通過できる何重かの障壁αの全部を通過する所要時間の確率密度を表わしている。動脈関数は、肺の障壁を通過して分散された関数と考えることができるため、このガンマ確率密度関数を、動脈カーブモデルとして適用し、動脈カーブモデルへの当てはめ（近似）を行うようにする。
【００５２】
動脈の時間濃度曲線ａ（t）は、上記の式（１１）で表されるガンマ確率密度関数の定数倍の関数でのカーブフィッティングによって、次式（１２）により表わされる。式（１２）において、ｔは、時刻を表し、t0，α，β，Ｋは、いずれもパラメータを表している。αは、障壁の数であり、βは、障壁αの全部を通過する所要時間であり、Ｋは、線形パラメータである。
［数１２］
a(t)=φ(t; t0, α, β, K) ・・・（１２）
【００５３】
動脈の時間濃度曲線算出部２ｃは、上記の式（１２）における４つのパラメータt0，α，β，Ｋの全てを変数とし、次式（１３）に示す、動脈カーブモデルと動脈の濃度値ａｉとの２乗誤差を最小にするt0，α，β，Ｋを算出する。実際には、t0には、造影剤が対象臓器に到達する前の特定の時刻を用い、α，β，Ｋを算出する方法を用いる場合が多い。
【数１３】
動脈の時間濃度曲線算出部２ｃは、カーブフィッティングした動脈の時間濃度曲線ａｉの算出値を、システム行列作成部２ｄに出力する。
【００５４】
システム行列作成部２ｄは、造影濃度曲線算出部２ｂから出力される組織の時間濃度曲線Ｃｉの値を入力するとともに、動脈の時間濃度曲線算出部２ｃから出力される動脈の時間濃度曲線ａｉの値を入力する。
【００５５】
システム行列作成部２ｄは、入力した組織の時間濃度曲線Ｃｉおよび動脈の時間濃度曲線ａｉから、解析対象となる先頭フェーズｉｂｓおよびフェーズ範囲ｉｂｓ≦ｉ≦Ｌ−１を決定する。解析対象となる先頭フェーズｉｂｓは、動脈および組織の造影濃度が上昇し始めるフェーズ、もしくは、それよりも前のフェーズを用いる。フェーズとは、解析上のサンプリング間隔により決まる時刻のことであり、フェーズ範囲とはいわゆる解析範囲のことである。
【００５６】
システム行列作成部２ｄは、決定したフェーズ範囲における動脈の時間濃度曲線ａｉを時間軸方向に離散化したもの、つまり、ｉｂｓ≦ｉ≦Ｌ−１のフェーズ範囲をＬ個の等区間に分割したＬ個のデータを、次式（１４）に示すようにして算出する。
［数１４］
Ａｉ＝ａ（ｉＴＲ） ・・・（１４）
【００５７】
上記の式（１４）において、ＴＲは、画像撮影の時間間隔（繰り返し時間）を表わしているが、これに限らず、ダイナミック時相、または、ダイナミック時相から得られた実際の解析時相でもよい。また、測定値ａｉをそのままＡｉとして使用するようにしてもよい。
【００５８】
システム行列作成部２ｄは、組織の時間濃度曲線Ｃｉを時間軸方向に離散化することにより得た、異なる時刻での組織の濃度値からなる集合の、時間軸方向の前方に、Ｈ個の要素を追加する。
【００５９】
これにより、先頭フェーズｉｂｓより前方にＨ個の要素が追加され、収集開始時刻または所定時刻（収集開始時刻からＮフェーズ目）を基準にして、異なる時刻での組織の濃度値を含む複数の要素からなる第１の集合が１次元（時間軸方向）に配列され、（Ｈ＋Ｌ）×１の列ベクトルが構成される。
【００６０】
このシステム行列作成部２ｄは、組織に対応する時間濃度曲線を、第１の時刻（収集開始時刻または収集開始時刻からＮフェーズ目）を基準にして離散化し、異なる時刻での組織における造影剤の濃度に対応する値を含む複数の要素を１次元に配列した第１の集合を構成する第１の構成手段として機能する。
【００６１】
また、システム行列作成部２ｄは、動脈の時間濃度曲線Ａｉを時間軸方向に離散化することにより得た、異なる時刻での動脈の濃度値からなる集合の、時間軸方向の後方に、Ｈ個の要素を追加する。これにより、異なる時刻での動脈の濃度値を含む複数の要素からなる第２の集合が１次元（時間軸方向）に配列され、（Ｌ＋Ｈ）×１の列ベクトルが構成される。そしてシステム行列作成部２ｄは、１次元に配列された動脈の濃度値の列ベクトルをさらに２次元（空間方向）に配列した（Ｌ＋Ｈ）×（Ｌ＋Ｈ）の係数行列を構成する。
【００６２】
これにより、解析対象の先頭フェーズｉｂｓを基準にして、異なる時刻での動脈の濃度値を含む複数の要素からなる第２の集合が２次元（時間軸方向と空間方向）に配列され、（Ｌ＋Ｈ）×（Ｌ＋Ｈ）の係数行列が構成される。
【００６３】
このシステム行列作成部２ｄは、動脈に対応する時間濃度曲線を、第１の時刻（収集開始時刻または収集開始時刻からＮフェーズ目）より後方の第２の時刻（動脈の時間濃度曲線が上昇し始める時刻）を基準にして離散化し、異なる時刻での動脈における造影剤の濃度に対応する値を含む複数の要素を２次元に配列した第２の集合を構成する第２の構成手段としても機能する。
【００６４】
システム行列作成部２ｄは、動脈の時間濃度曲線とインパルス応答関数のデコンボリューションが組織の時間濃度曲線であることから、以上のように構成した、異なる時刻での組織の濃度値を含む複数の要素からなる列ベクトルＣ、異なる時刻での動脈の濃度値を含む複数の要素からなる係数行列Ａ、および未知数Ｘ（インパルス応答関数）を用いて、次式（１５）に示すようなＣ＝ＡＸのシステム行列を作成する。式（１５）において、Ｃは、（Ｈ＋Ｌ）×１ベクトルであり、Ａは、（Ｌ＋Ｈ）×（Ｌ＋Ｈ）行列であり、Ｘは、（Ｌ＋Ｈ）×１ベクトルである。ｉｂｓは、解析対象となる先頭フェーズを表わしている。
【数１５】
【００６５】
上記の式（１５）において、Ｈは、Ｌ×Ｌ係数行列に追加する要素の個数を表わしている。例えば、撮影（収集）開始時刻から先頭フェーズまでの範囲をフェーズ間隔で分割した個数、または、それより多くても少なくても構わない。追加要素の個数Ｈの数を少なくすれば、解析対象となる行列のサイズが小さく処理時間が短縮されるが、少なすぎると時間方向のずれに対する依存性が増大することとなる。従って、撮影開始時刻から先頭フェーズまでの範囲をフェーズ間隔で分割した個数の値を用いれば、時間方向のずれに対する依存性を小さくするのに十分な値となる。
【００６６】
組織の時間濃度曲線を時間軸方向に離散化した複数の要素のうち、追加分の要素値Ｐｉ（ｉ＝０，１，・・・Ｈ−１）には、０の要素を入力（追加）するが、部分的あるいは全部をＣｉとして、測定された組織の時間濃度曲線の値を用いるようにしても良い。この場合、実際にはＬ個以上のフェーズの画像データが撮影されている必要がある。
【００６７】
動脈の時間濃度曲線を時間軸方向に離散化した要素値ａi,jは、次式（１６）により表わされる。式（１６）において、ｉ＝０，１，・・・Ｌ−１であり、ｊ＝０，１，...Ｌ−１である。
【数１６】
【００６８】
上記の式（１６）において、追加分の要素値ｑｉ，ｊには、０の要素を入力（追加）するが、部分的あるいは全部を次式（１７）に示すようにして、測定された動脈の時間濃度曲線の値を用いるようにしても良い。この場合、実際にはＬ個以上のフェーズの画像データが撮影されている必要がある。
【数１７】
システム行列作成部２ｄは、作成した上記の式（１５）をデコンボリューション部２ｅに出力する。
【００６９】
デコンボリューション部２ｅは、システム行列作成部２ｄから出力された係数行列を入力し、SVD法やフーリエ変換法などの方法を用いてデコンボリューションを行う逆畳み込み積分手段として機能する。デコンボリューション部２ｅは、デコンボリューションにより得られた未知数Ｘから、上記の式（７）に示すようにしてインパルス応答関数Ｒｉを算出する。
【００７０】
デコンボリューション部２ｅは、上記の式（７）で算出されたインパルス応答関数Ｒｉから、ピークフェーズ（インパルス応答関数の最大値）の存在範囲を検索し、検索したピークフェーズからピーク範囲（幅）を算出する。
図３は、ピークフェーズを説明する図である。
【００７１】
図３に示すように、第１の実施形態において算出されるインパルス応答関数Ｒｉのピークは、従来のstandard SVD法に比べ、ほぼＨ個分後方のフェーズとなる。従って、デコンボリューション部２ｅは、ピークフェーズが後方にずれていることを考慮しなければならないため、局所血流動態に関する情報の算出の前に、予めピーク範囲を決定する必要がある。
【００７２】
ピーク範囲の算出方法として、例えば、先頭フェーズから時間軸方向にピークフェーズを検索し、検索したピークフェーズの前後で負の値の極小値をとるフェーズまでをピーク範囲とする。なお、負の値の極小値ではなく、０の値をとるフェーズまでをピーク範囲とするようにしてもよい。
【００７３】
デコンボリューション部２ｅは、検索したピークフェーズおよびピーク範囲に基づいて、局所血流動態に関する情報である、平均通過時間MTT（Mean Transit Time）、局所脳血流量CBF（Cerebral Blood Flow）、および局所脳血液量CBV（Cerebral Blood Volume）を、次式（１８）に示すようにして算出する。つまり、インパルス応答曲線のピークの高さから局所脳血流量CBFが求まり、インパルス応答曲線のピーク範囲の面積から局所脳血液量CBVが求まり、インパル応答曲線のピーク範囲から平均通過時間MTTが求まる。
［数１８］
CBF＝MAX(Ri）
CBV＝SUM(Ri)
MTT＝SUM(Ri)/MAX(Ri)・・・（１８）
【００７４】
このデコンボリューション部２ｅは、組織の伝達関数（インパルス応答関数）に基づいて、血流動態に関する情報を算出する血流情報算出手段としても機能する。
【００７５】
またデコンボリューション部２ｅは、算出した平均通過時間MTT、局所脳血流量CBF、および局所脳血液量CBVを次式（１９）に示すようにして補正する。
［数１９］
CBV=(100h/ρ)TRSUM(Ri)
CBF(r)=60(100h/ρ)MAX(Ri)
MTT2(r)=TR{SUM(Ri)/MAX(Ri)} ・・・（１９）
【００７６】
上記の式（１９）において、ρは、脳の比重を表わし、ｈは、ヘマトクリット値（一定量の血液中に含まれる赤血球の割合）を表わしている。なお、局所脳血流量CBFは、次式（２０）の関係を用いて算出（補正）するようにしてもよい。
［数２０］
CBF＝CBV/MTT2 ・・・（２０）
【００７７】
次に、図４のフローチャートを参照して、第１の実施形態における組織の血流動態に関する情報の算出処理について説明する。この処理では、特に、脳組織を例に挙げ、その血流動態に関する情報を算出する処理について説明する。
ステップＳ１において、時系列データ取得部２ａは、医用画像撮影装置が所定時間撮影して生成した時系列の撮影データを取得する。
【００７８】
ステップＳ２において、造影濃度曲線算出部２ｂは、ステップＳ１の処理で取得された時系列の撮影データから、各画素の画素値の変化を計測して造影濃度曲線を算出する。この造影濃度曲線を、組織の時間濃度曲線Ｃｉとする。
【００７９】
ステップＳ３において、動脈の時間濃度曲線算出部２ｃは、ステップＳ２の処理で算出された組織の時間濃度曲線Ｃｉから、CT画像やMRI画像の解析単位となる動脈領域を設定し、その動脈領域での時間濃度曲線ａiを算出する。このとき、動脈の時間濃度曲線算出部２ｃは、取得した撮影データがMRI画像の場合、ガンマ確率密度関数の定数倍の関数を動脈カーブモデルとして適用し、算出した動脈の時間濃度曲線ａｉのカーブフィッティングを行う。
ステップＳ４において、システム行列作成部２ｄは、解析対象となる先頭フェーズｉｂｓおよびフェーズ範囲ｉｂｓ≦ｉ≦Ｌ−１を決定する。
【００８０】
ステップＳ５において、システム行列作成部２ｄは、ステップＳ２の処理で算出された組織の時間濃度曲線Ｃｉのうち、ステップＳ４の処理で決定したフェーズ範囲ｉｂｓ≦ｉ≦Ｌ−１における組織の時間濃度曲線を時間軸方向に離散化し、離散化することにより得た異なる時刻での組織の濃度値からなる集合を１次元に配列して列ベクトルを構成する。
【００８１】
また、ステップＳ５において、システム行列作成部２ｄは、ステップＳ３の処理で算出された動脈の時間濃度曲線ａiのうち、ステップＳ４の処理で決定したフェーズ範囲ｉｂｓ≦ｉ≦Ｌ−１における動脈の時間濃度曲線を時間軸方向に離散化し、離散化することにより得た異なる時刻での動脈の濃度値からなる集合を２次元に配列して係数行列を構成する。
【００８２】
ステップＳ６において、システム行列作成部２ｄは、ステップＳ５の処理で構成された異なる時刻での組織の濃度値からなる列ベクトル（集合）の、時間軸方向の前方にＨ個の要素を追加する。これにより、解析対象の先頭フェーズｉｂｓが時間軸方向にＨ個分ずらされ、収集開始時刻または所定時刻を基準にして、異なる時刻での組織の濃度値を含む複数の要素からなる第１の集合が１次元に配列され、（Ｈ＋Ｌ）×１の列ベクトルＣが構成される。
【００８３】
また、ステップＳ６において、システム行列作成部２ｄは、ステップＳ５の処理で構成された異なる時刻での動脈の濃度値からなる係数行列（集合）の、時間軸方向の後方にＨ個の要素を追加する。これにより、動脈の造影濃度が上昇し始めるフェーズもしくは、それよりも前のフェーズを基準にして、異なる時刻での動脈の濃度値を含む複数の要素からなる第２の集合が２次元に配列され、（Ｌ＋Ｈ）×（Ｌ＋Ｈ）の係数行列Ａが構成される。
【００８４】
ステップＳ７において、システム行列作成部２ｄは、ステップＳ６の処理で得られた、異なる時刻での組織の濃度値を含む複数の要素からなる列ベクトルＣ、異なる時刻での動脈の濃度値を含む複数の要素からなる係数行列Ａ、および、未知数Ｘ（インパルス応答関数）を用いて、上記の式（１５）に示すシステム行列を作成する。
図５は、第１の実施形態に係る式（１５）のシステム行列を模式的に示す図である。
【００８５】
図５において、Ｃは、組織の時間濃度曲線Ｃｉを時間軸方向に離散化した（Ｈ＋Ｌ）×１の列ベクトルを模式的に表わしたものであり、Ａは、動脈の時間濃度曲線ａｉを時間軸方向に離散化した（Ｌ＋Ｈ）×（Ｌ＋Ｈ）の係数行列を模式的に表わしたものであり、Ｘは未知数（インパルス応答関数）を模式的に表わしたものである。
【００８６】
同図に示すように、組織の時間濃度曲線を構成する列ベクトルＣにおいて、前方にＨ個の要素ｐｉが追加され、動脈の時間濃度曲線を構成する係数行列Ａにおいて、後方にＨ個の要素ｑｉ，ｊが追加されている。追加要素の個数Ｈには、撮影開始時刻から先頭フェーズｉｂｓまでの範囲をフェーズ間隔で分割した個数の値を用いることが好ましい。追加分の要素値Ｐｉ（ｉ＝０，１，・・・Ｈ−１）と要素値ｑｉ，ｊには、０の要素が入力（追加）されるが、部分的あるいは全部に、実際の測定値を用いるようにしても良い。
【００８７】
このようなシステム行列を構成することで、従来のblock circulant SVD法のように係数行列を２倍にする必要がなくなり、デコンボリューションの高速処理化が可能になる。
【００８８】
図４の説明に戻る。ステップＳ８において、デコンボリューション部２ｅは、ステップＳ７の処理で作成されたシステム行列を、SVD法あるいはフーリエ変換法などの方法を用いてデコンボリューションを行い、ベクトルＸを求め、上記の式（７）に示すようにしてインパルス応答関数Ｒｉを算出する。
【００８９】
ステップＳ９において、デコンボリューション部２ｅは、ステップＳ８の処理で算出されたインパルス応答関数Ｒｉから、ピークフェーズを検索し、検索したピークフェーズからピーク範囲を算出する。
【００９０】
ステップＳ１０において、デコンボリューション部２ｅは、ステップＳ９の処理で算出されたピークフェーズおよびピーク範囲に基づいて、局所血流動態に関する情報である、平均通過時間MTT、局所脳血流量CBF、および局所脳血液量CBVを上記の式（１８）に従って算出する。
【００９１】
ステップＳ１１において、デコンボリューション部２ｅは、ステップＳ１０の処理で算出された平均通過時間MTT、局所脳血流量CBF、および局所脳血液量CBVを、上記の式（１９）に示すようにして補正する。
【００９２】
以上のように、血流動態解析装置１は、解析対象となる組織の時間濃度曲線を構成する列ベクトルの前方に所定個数の要素を追加し、動脈の時間濃度曲線を構成する係数行列の後方に所定個数の要素を追加し、組織の時間濃度曲線の先頭フェーズと、動脈の時間濃度曲線の先頭フェーズの基準位置をずらすことによって、処理時間を増大させることなく、かつ、時間方向のずれに対する依存性を小さくして、血管の血流動態に関する情報である、組織血流量、組織血液量、および平均通過時間を安定に算出することが可能となる。
【００９３】
これにより、病院側は、多数の患者の検査を実施することができるようになり、患者の検査待ち時間を短縮することができる。また、安定に解析結果を得ることができるため、再検査（再解析）といった事態を減らすことができ、患者への負担を軽減することができる。
【００９４】
また、血流動態解析装置１は、算出した平均通過時間MTT、局所脳血液量CBV、および局所脳血流量CBFに基づいて、対応する画像を表示部１ｈに表示させることができ、操作者は、血流動態の状態を容易に把握することが可能となる。
【００９５】
以上においては、対象組織として、脳組織を例に挙げ説明したが、第１の実施形態では、脳組織に限られるものではなく、肝臓、心臓、あるいは肺などの他の組織を対象とすることも勿論可能である。
【００９６】
また以上において、脳組織についてヘマトクリット値の補正および脳の比重による補正を行うようにした。そこで、例えば、肝臓や肺などの呼吸による動きがある組織、あるいは、心臓を対象組織とした場合については、組織の動き補正を行う必要がある。また動き補正だけでなく、ヘマトクリット値の補正、その他の補正を適宜行う必要がある。
【００９７】
（第２の実施形態）
図６は、本発明の第２の実施形態に係る流体解析装置の構成例を示す機能ブロック図である。
【００９８】
流体解析装置３は、例えば、MRI装置４に内蔵される。但し、ネットワークを介して流体解析装置３をMRI装置４或いはMRI装置４において収集された画像データを保管する画像処理装置や画像サーバに接続してもよい。
【００９９】
MRI装置４は、イメージング部４ａを備える。イメージング部４ａは、被検体の流体の非造影イメージング又は造影イメージングを実行することによって、複数フレーム分の非造影流体MR画像データ又は造影流体MR画像データを収集する機能を有する。流体としては、血管内の血流や脳脊髄液(CSF: cerebrospinal fluid)が挙げられる。MRI装置を用いた血管の撮影は、磁気共鳴血管撮影法(MRA: magnetic resonance angiography)と呼ばれる。
【０１００】
具体的には、イメージング部４ａは、静磁場磁石、傾斜磁場コイル、傾斜磁場電源、送信用高周波(RF: radio frequency)コイル、受信用RFコイル、送信器、受信器、シーケンサ及びコンピュータ等の構成要素を有する。そして、イメージング部４ａでは、コンピュータにおいて設定されたパルスシーケンス等の撮像条件に従って、静磁場磁石内の静磁場下にセットされた被検体に傾斜磁場コイルからは傾斜磁場が、送信用コイルからはRFパルスが、それぞれ印加される。そして、被検体から発生した核磁気共鳴(NMR: nuclear magnetic resonance)信号が受信用コイルで受信され、受信器においてデジタル化されたNMR信号に対する画像再構成処理がコンピュータにおいて実行されることによってMR画像データが生成される。また、傾斜磁場コイル及び送信用RFコイルには、それぞれコンピュータからシーケンサを通じて出力されるパルスシーケンスに従って傾斜磁場電源及び送信器から傾斜磁場パルス及びRFパルスが印加される。これにより、傾斜磁場コイル及び送信用RFコイルからは、所望の波形の傾斜磁場パルス及びRFパルスが被検体に印加される。
図７は、図６に示すMRI装置４のイメージング部４ａにおいて実行される非造影流体イメージング用の撮像条件の第１の例を示す図である。
【０１０１】
図７において横軸は時間を示す。図７は、非造影で周期性のないCSFを描出するためのパルスシーケンスの一例を示している。図７に示すように、ラベリングパルス（タグ付けパルスとも言う）をCSFまたは着目領域に印加し、続いて連続的に複数回に亘ってイメージングデータの収集(DATA ACQUISITION 1, DATA ACQUISITION 2, DATA ACQUISITION 3,...)を繰り返すダイナミックスキャンによって、異なる時刻t1, t2, t3, ...に対応する複数フレーム分の時系列の非造影流体画像データを収集することができる。
【０１０２】
ラベリングパルスとしては、t-SLIP (Time-SLIP: Time Spatial Labeling Inversion Pulse)、飽和(SAT: saturation)パルス（Rest grid pulseとも呼ばれる）、SPAMM (spatial modulation of magnetization)パルス及びダンテ(DANTE)パルスが知られている。同一種類のラベリングパルスを複数回印加したり、異種の複数のラベリングパルスを組み合わせて印加してもよい。
【０１０３】
例えば、複数の領域選択的90°SATパルスを組み合わせると、選択スラブ領域に放射状やストライプ状のパターンで飽和された領域を形成することができる。また、SPAMMパルスやダンテパルスの印加によってもストライプパターン、グリッドパターン（格子状のパターン）、放射状のパターン等の所望のパターンで飽和された領域を形成することができる。このため、パターンがCSFの流れに伴って移動する複数フレーム分のダイナミックCSF画像データを収集することができる。
【０１０４】
また、t-SLIPを用いたイメージングには、flow-in法とflow-out法がある。flow-in法は、t-SLIPを構成する領域選択反転回復(IR: inversion recovery)パルスを着目領域に印加し、着目領域の外部から着目領域に流入するラベリングされていないCSFを描出する方法である。一方、flow-out法は、t-SLIPを構成する領域非選択IRパルスを印加するとともにラベリング領域に領域選択IRパルスを印加し、ラベリング領域から着目領域に流入するラベリングされたCSFを描出する方法である。
【０１０５】
ラベリングパルスは、ECG (electro cardiogram)信号、呼吸センサによる同期信号又は脈波同期(PPG: peripheral pulse gating)信号等の生体信号に同期させて印加してもよい。
【０１０６】
一方、イメージングデータの収集用のシーケンスとしては、SSFP (steady state free precession)シーケンスやFASE (FastASE: fast asymmetric spin echoまたはfast advanced spin echo)シーケンス等の任意のシーケンスが使用できる。
図８は、図６に示すMRI装置４のイメージング部４ａにおいて実行される非造影流体イメージング用の撮像条件の第２の例を示す図である。
【０１０７】
図８において横軸は時間を示し、縦方向も時間を表している。図８は、IRパルスを用いて非造影で被検体の血流やCSF等の流体をイメージングするためのパルスシーケンスの一例を示している。すなわち、IRパルスの印加時刻からイメージングデータの収集タイミングまでの反転時間(TI: inversion time)をTI1, TI2, TI3, ...と徐々に変えながら、IRパルスの印加及びイメージングデータの収集(DATA ACQUISITION 1, DATA ACQUISITION 2, DATA ACQUISITION 3,...)が繰り返し実行される。
【０１０８】
IRパルスは単一のパルスに限らず、同時とみなせる十分に短い間隔で異種又は同種の複数のIRパルスを印加してもよい。領域選択的IRパルスは、上述のflow-in法の場合には撮像領域に印加され、flow-out法の場合には撮像領域の上流領域に印加される。そして、領域選択的IRパルス、周波数選択的IRパルス、非領域選択的IRパルス及び非周波数選択的IRパルスの組み合わせによって所望の対象物の磁化を倒すことによって撮像領域又は流体がラベリングされる。そして、TIに応じた距離だけ流体が移動したタイミングで撮像領域から収集したイメージングデータの画像再構成処理によって、撮像領域に流入する流体が選択的に描出された流体画像データを生成することができる。
【０１０９】
これにより、異なる複数のTIに対応する複数フレーム分の非造影流体画像データを収集することができる。すなわち、複数のTIに対応する複数フレーム分の画像データをTIの短い順に並べると、ダイナミック画像データのようにTIの増加に伴って流体が徐々に流れていく様子が描出される連続的な画像データを生成することができる。
【０１１０】
尚、IRパルスの印加及びイメージングデータの収集は、データ収集タイミングにおける流体の速度をより一定に近づけるために、ECG信号等の生体信号に同期させてもよい。生体信号を用いた同期イメージングは、血流のように周期性のある流体のイメージングに好適である。
【０１１１】
また、イメージングデータを収集するためのパルスシーケンスとしては、SE (spin echo)シーケンス、FSE (fast spin echo)シーケンス、FASEシーケンス、EPI (echo planar imaging) シーケンス等のパルスシーケンスが用いられる。
図９は、図６に示すMRI装置４のイメージング部４ａにおいて実行される非造影流体イメージング用の撮像条件の第３の例を示す図である。
【０１１２】
領域選択的IRパルスの印加及び図９に示すような血流を含む実線で示すイメージング領域からのイメージングデータの収集を繰り返し、かつTIを一定にしつつ領域選択的IRパルスの印加領域を点線で示すようにTAGGED REGION 1,TAGGED REGION 2, TAGGED REGION 3, ...に徐々に変えるイメージングによってもダイナミック画像データのように流体が徐々に流れていく様子が描出される連続的な画像データを生成することができる。すなわち、領域選択的IRパルスの印加領域に存在していた血流はラベリングされ、領域選択的IRパルスの印加領域から下流に向かってTIに応じた距離だけ流出したタイミングでイメージングデータが収集される。従って、ラベリングされた血流が選択的に描出された画像データを生成することができる。
【０１１３】
このため、領域選択的IRパルスの印加領域を変えて繰り返しイメージングデータを収集することによって、異なる複数のラベリング領域に対応し、ラベリングされた血流がイメージング領域の異なる位置に存在する複数フレーム分の非造影流体画像データを収集することができる。
【０１１４】
さらに、血管のより上流において血流が描出された画像データを各ラベリング領域に対応する画像データと最大値投影(MIP: maximum intensity projection)処理等の合成処理によって合成することにより、血流が徐々に流れていく様子が描出される連続的な画像データを生成することができる。
尚、図８及び図９においてIRパルスの代わりにSATパルスを用いてもよい。
【０１１５】
ここまでは、非造影イメージングの撮像条件について例を示したが、第１の実施形態と同様に公知の造影イメージング法によって造影画像データを収集することもできる。この場合、複数の時相に対応する時系列の複数フレーム分の造影画像データが収集される。
【０１１６】
一方、流体解析装置３は、イメージング部４ａにおいて収集された非造影流体MR画像データ又は造影流体MR画像データに基づいて、流体の動態を表す情報を取得する機能を有する。流体解析装置３は、例えば、コンピュータの演算装置に流体解析プログラムを読み込ませることにより、コンピュータをイメージデータ取得部３ａ、組織画素変化曲線算出部３ｂ、流体画素変化曲線算出部３ｃ、システム行列作成部３ｄ及びデコンボリューション部３ｅとして機能させたものである。但し、回路を用いて流体解析装置３の一部又は全部を構成してもよい。
【０１１７】
イメージデータ取得部３ａは、イメージング部４ａにおいて収集された複数フレーム分の連続的な流体画像データを取得して組織画素変化曲線算出部３ｂに与える機能を有する。
【０１１８】
組織画素変化曲線算出部３ｂは、複数フレーム分の連続的な流体画像データにおいて、組織に対応する複数の画素におけるフレーム間における画素値の変化から、組織内を通過するラベリングされた血流等の流体の濃度に応じた量の変化を求め、それぞれ組織の画素濃度変化曲線として算出する機能と、算出した複数の組織の画素濃度変化曲線を流体画素変化曲線算出部３ｃ及びシステム行列作成部３ｄに出力する機能を有する。
【０１１９】
尚、非造影イメージング法では、実際には造影剤が注入されない。このため、非造影イメージング法の場合には、血流等の流体の成分の量に関する濃度ではなく、ラベリングされた流体成分の流体全量に対する比に対応する論理的な濃度が求められる。
【０１２０】
例えば、時刻(t1, t2, t3, ...)、TI (TI1, TI2, TI3, ...)或いはラベリング領域(TAGGED REGION 1,TAGGED REGION 2, TAGGED REGION 3, ...)等のフレーム間で変えられるパラメータのi番目の値に対応する画素値をSiとし、初期値(t1, TI1, TAGGED REGION 1)等の基準となるパラメータの値に対応する画素値をS0とする。そうすると、第１の実施形態で示した式(１０)によって、組織における画素濃度に比例する量の変化を画素濃度変化曲線として算出することができる。或いは、フレーム間で変えられるパラメータに対応する組織の画素値S自体の変化を表す時間強度曲線(TIC: time intensity curve)を画素濃度変化曲線としてもよい。
画素濃度に比例する量及び画素値は、流体の流れに対応する量となる。従って、組織における画素濃度変化曲線は、組織における流体の動態情報を含んでいる。
【０１２１】
流体画素変化曲線算出部３ｃは、組織の複数の画素における画素濃度変化曲線に基づいて、動脈等の流体領域内を循環する流体の画素濃度変化曲線を算出する機能と、算出した流体の画素濃度変化曲線をシステム行列作成部３ｄに出力する機能を有する。流体の画素濃度変化曲線は、流体領域内を循環する流体に対応する画素値の濃度に応じた量のフレーム間の変化を表す曲線である。
【０１２２】
流体の画素濃度変化曲線は、第１の実施形態と同様な方法で算出することができる。このため、第１の実施形態と同様に組織の画素濃度曲線は、複数の画素についてそれぞれ求められるのに対して、流体の画素濃度曲線は流体の画素値の濃度に応じた量を代表する値を用いて１つの曲線として求められる。従って、流体部位として複数の画素を含む関心領域(ROI: region of interest)が設定された場合には、第１の実施形態と同様に画素値の濃度に応じた量を代表する値が算出される。代表値としては、画素間の平均値を用いる方法や式（１１）に示すようなフィッティングパラメータを有するガンマ確率密度関数の定数倍の関数を用いたカーブフィッティングにより算出された値を用いる方法がある。流体領域に対応する画素についても、第１の実施形態と同様に自動的あるいはユーザからの指定によって特定することができる。
【０１２３】
システム行列作成部３ｄは、組織の各画素に対応する画素濃度変化曲線に基づいて、時刻(t1, t2, t3, ...)、TI (TI1, TI2, TI3, ...)或いはラベリング領域(TAGGED REGION 1,TAGGED REGION 2, TAGGED REGION 3, ...)等のパラメータの値ごとの組織の画素値の濃度に応じた量を要素とする１次元の列ベクトルＣを作成する機能と、流体の画素濃度変化曲線に基づいて、パラメータの値ごとの流体の画素値の濃度に応じた量を要素とする１次元の列ベクトルを作成する機能とを有する。加えて、システム行列作成部３ｄは、組織に対応する列ベクトルＣ、流体に対応する列ベクトルを２次元に配列した係数行列Ａ及び未知数Ｘを用いて、Ｃ＝ＡＸのシステム行列を作成する機能と、作成したシステム行列をデコンボリューション部３ｅに出力する機能とを有する。
【０１２４】
組織に対応する列ベクトルＣ、流体に対応する列ベクトル、係数行列Ａ及びシステム行列Ｃ＝ＡＸは、それぞれ第１の実施形態における方法と同様な方法で作成することができる。換言すれば以下のような説明になる。
図１０は、図６に示すシステム行列作成部３ｄにおける解析対象となる要素の作成方法を説明する図である。
【０１２５】
図１０において横軸は、時刻(t1, t2, t3, ...)、TI (TI1, TI2, TI3, ...)或いはラベリング領域(TAGGED REGION 1,TAGGED REGION 2, TAGGED REGION 3, ...)等のパラメータ値を示し、縦軸は画素値の濃度を示す。
【０１２６】
図１０に示すように、組織の各画素に対応する画素濃度曲線C1i, C2i, C3i, ...及び流体の画素濃度曲線aiが得られている場合に、複数の組織の画素濃度変化曲線C1i, C2i, C3i, ...がそれぞれ少なくとも立ち上がる前となる第１の点及び流体の画素濃度変化曲線aiが少なくとも立ち上がる前となる第２の点が決定される。例えば、第１の点は、時刻(t1, t2, t3, ...)、TI (TI1, TI2, TI3, ...)或いはラベリング領域(TAGGED REGION 1,TAGGED REGION 2, TAGGED REGION 3, ...)等のイメージングにおいて変えられたパラメータの初期値(t1, TI1, TAGGED REGION 1)に決定される。また、第２の点は、第１の点よりも順番が後の値に決定される。
【０１２７】
さらに、組織の画素濃度曲線C1i, C2i, C3i, ...に対応する解析対象となる要素と流体の画素濃度曲線aiに対応する解析対象となる要素が同じ数になるように、第１の点及び第２の点を開始点とする同じ広さの解析対象範囲RANGE1, RANGE2が組織の画素濃度曲線C1i, C2i, C3i, ...及び流体の画素濃度曲線aiに対してそれぞれ設定される。尚、組織の画素濃度曲線C1i, C2i, C3i, ...の解析対象範囲RANGE1のみに含まれるパラメータ値における前方の一部又は全部の要素をゼロとしてもよい。一方、流体の画素濃度曲線aiの解析対象範囲RANGE2のみに含まれるパラメータ値における後方の一部又は全部の要素をゼロとしてもよい。
【０１２８】
この結果、第１の点における組織の画素値の濃度に応じた量又はゼロが組織の画素濃度曲線C1i, C2i, C3i, ...に対応する解析対象となる最初の要素となり、第２の点における流体の画素値の濃度に応じた量が流体の画素濃度曲線aiに対応する解析対象となる最初の要素となる。
【０１２９】
このように作成された同じ数の組織の画素濃度曲線に対応する要素及び流体の画素濃度曲線に対応する要素を用いて式（１５）に示すようなシステム行列を作成することができる。
【０１３０】
デコンボリューション部３ｅは、第１の実施形態におけるデコンボリューション部２ｅが有する機能と同様の機能を有する。すなわち、デコンボリューション部３ｅは、組織に対応する列ベクトルＣと流体に対応する列ベクトルを２次元に配列した係数行列Ａの逆畳み込み積分によって未知数Ｘを算出する機能と、未知数Ｘから得られるインパルス応答関数に基づいて流体の平均通過時間MTT、局所脳血流量CBF、局所脳血液量CBV、局所CSF流量、局所CSF量等の流体動態に関する情報を算出する機能を有する。
【０１３１】
また、必要に応じてデコンボリューション部３ｅは、第１の実施形態におけるデコンボリューション部２ｅと同様に、流体動態に関する情報を式（１９）に示すように補正するように構成される。
次に、流体解析装置３の動作及び作用について説明する。
【０１３２】
図１１は、図６に示す流体解析装置３におけるデータ処理の流れを示すフローチャートである。尚、第１の実施形態において図４に示されている詳細な処理と同様な処理については省略する。
まず予め、MRI装置４のイメージング部４ａは非造影イメージング又は造影イメージングを実行し、複数フレーム分のMR流体画像データを収集する。
【０１３３】
そして、ステップＳ２１において、流体解析装置３のイメージデータ取得部３ａは、イメージング部４ａから複数フレーム分の流体画像データを取得して組織画素変化曲線算出部３ｂに与える。
【０１３４】
次に、ステップＳ２２において、組織画素変化曲線算出部３ｂは、複数フレーム分の流体画像データの組織に対応する複数の画素におけるフレーム間における画素値の濃度に応じた量の変化を求めることによって組織の画素濃度変化曲線を算出する。
【０１３５】
次に、ステップＳ２３において、流体画素変化曲線算出部３ｃは、組織の複数の画素における画素濃度変化曲線に基づいて、流体の画素濃度変化曲線を算出する。
【０１３６】
次に、ステップＳ２４において、システム行列作成部３ｄは、組織の画素濃度変化曲線及び流体の画素濃度変化曲線に基づいて、組織の画素濃度に対応する列ベクトルをＣ、流体の画素濃度に対応する列ベクトルを２次元に配列した係数行列をＡ、未知数をＸとするＣ＝ＡＸのシステム行列を作成する。このとき、組織に対応する列ベクトルＣの最初の要素は、組織の画素濃度変化曲線の横軸上の第１のパラメータ値に対応する値とされる。一方、流体に対応する列ベクトルの最初の要素は、流体の画素濃度変化曲線の横軸上の第１のパラメータ値よりも後方の第２のパラメータ値に対応する値とされる。また、組織に対応する列ベクトルＣの要素数と流体に対応する列ベクトルの要素数が同じになるように各列ベクトルが作成される。そのために、列ベクトルを構成する要素として画素濃度の代わりにゼロを用いてもよい。
【０１３７】
次に、ステップＳ２５において、デコンボリューション部３ｅは、システム行列の未知数Ｘを算出し、未知数Ｘから得られるインパルス応答関数に基づいて流体動態に関する情報を算出する。この流体動態に関する情報は、式（１９）によって補正してもよい。
【０１３８】
このように、流体解析装置３では、複数の時相に対応する時系列の造影画像データに限らず、撮像パラメータを変えて収集された連続的な複数フレーム分の非造影画像データに基づいて、造影画像データに対する処理と同様な処理で血流やCSF等の流体の動態に関する情報を求めることができる。
【０１３９】
尚、第２の実施形態では、流体解析装置３がMRI装置で収集された画像データに基づいて流体動態に関する情報を取得する例について説明したが、Ｘ線CT装置、超音波診断装置又はPET装置等の医用画像撮影装置で収集された造影画像データ又は非造影画像データに基づいて流体解析装置３が流体動態に関する情報を取得するように構成することもできる。
【０１４０】
（他の実施形態）
以上、特定の実施形態について記載したが、記載された実施形態は一例に過ぎず、発明の範囲を限定するものではない。ここに記載された新規な方法及び装置は、様々な他の様式で具現化することができる。また、ここに記載された方法及び装置の様式において、発明の要旨から逸脱しない範囲で、種々の省略、置換及び変更を行うことができる。添付された請求の範囲及びその均等物は、発明の範囲及び要旨に包含されているものとして、そのような種々の様式及び変形例を含んでいる。
【符号の説明】
【０１４１】
１ 血流動態解析装置
２ 局所血流解析アプリケーション
２ａ 局所血流解析アプリケーション
２ｂ 造影濃度曲線算出部
２ｃ 動脈の時間濃度曲線算出部
２ｄ システム行列作成部
２ｅ デコンボリューション部
３ 流体解析装置
３ａ イメージデータ取得部
３ｂ 組織画素変化曲線算出部
３ｃ 流体画素変化曲線算出部
３ｄ システム行列作成部
３ｅ デコンボリューション部
４ MRI装置
４ａ イメージング部
【特許請求の範囲】
【請求項１】
医用画像撮影装置により得られた複数の時相における被検体の造影画像データに基づいて、前記複数の時相間における各画素の画素値の変化を計測して組織における造影剤の濃度に対応する量の時間変化を表す時間濃度曲線を算出する組織の時間濃度曲線算出手段と、
前記組織に対応する前記時間濃度曲線から動脈領域を設定し、動脈に対応する時間濃度曲線を算出する動脈の時間濃度曲線算出手段と、
前記組織に対応する前記時間濃度曲線を、第１の時刻を基準にして離散化し、異なる時刻での前記組織における前記造影剤の濃度に対応する値を含む複数の要素を１次元に配列した第１の集合を構成する第１の構成手段と、
前記動脈に対応する前記時間濃度曲線を、前記第１の時刻より後方の第２の時刻を基準にして離散化し、異なる時刻での前記動脈における前記造影剤の濃度に対応する値を含む複数の要素を２次元に配列した第２の集合を構成する第２の構成手段と、
前記第１の集合と前記第２の集合の逆畳み込み積分を行い、前記組織の伝達関数を算出する逆畳み込み積分手段と、
前記伝達関数に基づいて、血流動態に関する情報を算出する血流情報算出手段と
を備える血流動態解析装置。
【請求項２】
前記第１の時刻は、前記画像データの収集開始時刻、または前記収集開始時刻と前記第２の時刻の間の所定時刻であり、前記第２の時刻は、前記動脈に対応する前記時間濃度曲線が上昇し始める時刻である請求項１に記載の血流動態解析装置。
【請求項３】
前記第１の集合は、１次元に配列されたＬ個の異なる時刻での前記組織における前記造影剤の濃度に対応する値の集合に対して、時間軸方向の前方にＨ個の第１の要素が追加された集合であり、前記第２の集合は、２次元に配列されたＬ個の異なる時刻での前記動脈における前記造影剤の濃度に対応する値の集合に対して、前記時間軸方向の後方にＨ個の第２の要素が追加された集合である請求項１又は２に記載の血流動態解析装置。
【請求項４】
前記第１の要素は、０及び前記組織における前記造影剤の濃度に対応する値の少なくとも一方であり、前記第２の要素は、０及び前記動脈における前記造影剤の濃度に対応する値の少なくとも一方である請求項３に記載の血流動態解析装置。
【請求項５】
前記血流情報算出手段は、前記伝達関数の最大値を検索し、検索した前記伝達関数の最大値からピーク範囲を算出し、前記伝達関数の最大値および前記ピーク範囲に基づいて、前記血流動態に関する情報を算出する請求項１乃至４のいずれか１項に記載の血流動態解析装置。
【請求項６】
コンピュータを、
医用画像撮影装置により得られた複数の時相における被検体の造影画像データに基づいて、前記複数の時相間における各画素の画素値の変化を計測して組織における造影剤の濃度に対応する量の時間変化を表す時間濃度曲線を算出する組織の時間濃度曲線算出手段、
前記組織に対応する前記時間濃度曲線から動脈領域を設定し、動脈に対応する時間濃度曲線を算出する動脈の時間濃度曲線算出手段、
前記組織に対応する前記時間濃度曲線を、第１の時刻を基準にして離散化し、異なる時刻での前記組織における前記造影剤の濃度に対応する値を含む複数の要素を１次元に配列した第１の集合を構成する第１の構成手段、
前記動脈に対応する前記時間濃度曲線を、前記第１の時刻より後方の第２の時刻を基準にして離散化し、異なる時刻での前記動脈における前記造影剤の濃度に対応する値を含む複数の要素を２次元に配列した第２の集合を構成する第２の構成手段、
前記第１の集合と前記第２の集合の逆畳み込み積分を行い、前記組織の伝達関数を算出する逆畳み込み積分手段、及び
前記伝達関数に基づいて、血流動態に関する情報を算出する血流情報算出手段、
として機能させる血流動態解析プログラム。
【請求項７】
被検体の複数のフレームの画像データに基づいて、組織に対応する画素において画素値の濃度に応じた量の前記複数のフレーム間における第１の変化を算出する一方、流体に対応する画素値の濃度に応じた量の前記複数のフレーム間における第２の変化を算出する濃度変化算出手段と、
前記第１の変化に基づいて前記第１の変化の第１の点に対応する要素を最初の要素とする解析対象となる数の第１の複数の要素を作成する一方、前記第２の変化に基づいて前記第２の変化の前記第１の点よりも後の第２の点に対応する要素を最初の要素とする前記解析対象となる数の第２の複数の要素を作成し、前記第１の複数の要素及び前記第２の複数の要素に基づいて逆畳み込み積分処理を伴って流体の動態に関する情報を算出する流体情報算出手段と、
を備える流体解析装置。
【請求項８】
前記濃度変化算出手段は、複数のフレームの造影画像データに基づいて、造影剤の濃度に応じた量の時間変化を前記第１及び第２の変化として算出するように構成される請求項７に記載の流体解析装置。
【請求項９】
前記濃度変化算出手段は、複数のフレームの非造影画像データに基づいて、前記第１及び第２の変化を算出するように構成される請求項７に記載の流体解析装置。
【請求項１０】
前記濃度変化算出手段は、磁気共鳴イメージング装置により共通のラベリングパルスからイメージングデータの収集タイミングまでの時間を変えてダイナミック収集された複数のフレームの非造影画像データに基づいて、前記第１及び第２の変化を算出するように構成される請求項９に記載の流体解析装置。
【請求項１１】
前記濃度変化算出手段は、磁気共鳴イメージング装置により反転回復パルス又は飽和パルスの印加とイメージングデータの収集とを繰り返し、かつ前記反転回復パルス又は飽和パルスの印加タイミングから前記イメージングデータの収集タイミングまでの時間を変えて収集された複数のフレームの非造影画像データに基づいて、前記第１及び第２の変化を算出するように構成される請求項９に記載の流体解析装置。
【請求項１２】
前記濃度変化算出手段は、磁気共鳴イメージング装置により反転回復パルス又は飽和パルスの印加とイメージングデータの収集とを繰り返し、かつ前記反転回復パルス又は飽和パルスの印加領域を変えて収集された複数のフレームの非造影画像データに基づいて、前記第１及び第２の変化を算出するように構成される請求項９に記載の流体解析装置。
【請求項１３】
コンピュータを、
被検体の複数のフレームの画像データに基づいて、組織に対応する画素において画素値の濃度に応じた量の前記複数のフレーム間における第１の変化を算出する一方、流体に対応する画素値の濃度に応じた量の前記複数のフレーム間における第２の変化を算出する濃度変化算出手段、及び
前記第１の変化に基づいて前記第１の変化の第１の点に対応する要素を最初の要素とする解析対象となる数の第１の複数の要素を作成する一方、前記第２の変化に基づいて前記第２の変化の前記第１の点よりも後の第２の点に対応する要素を最初の要素とする前記解析対象となる数の第２の複数の要素を作成し、前記第１の複数の要素及び前記第２の複数の要素に基づいて逆畳み込み積分処理を伴って流体の動態に関する情報を算出する流体情報算出手段、
として機能させる流体解析プログラム。
【請求項１】
医用画像撮影装置により得られた複数の時相における被検体の造影画像データに基づいて、前記複数の時相間における各画素の画素値の変化を計測して組織における造影剤の濃度に対応する量の時間変化を表す時間濃度曲線を算出する組織の時間濃度曲線算出手段と、
前記組織に対応する前記時間濃度曲線から動脈領域を設定し、動脈に対応する時間濃度曲線を算出する動脈の時間濃度曲線算出手段と、
前記組織に対応する前記時間濃度曲線を、第１の時刻を基準にして離散化し、異なる時刻での前記組織における前記造影剤の濃度に対応する値を含む複数の要素を１次元に配列した第１の集合を構成する第１の構成手段と、
前記動脈に対応する前記時間濃度曲線を、前記第１の時刻より後方の第２の時刻を基準にして離散化し、異なる時刻での前記動脈における前記造影剤の濃度に対応する値を含む複数の要素を２次元に配列した第２の集合を構成する第２の構成手段と、
前記第１の集合と前記第２の集合の逆畳み込み積分を行い、前記組織の伝達関数を算出する逆畳み込み積分手段と、
前記伝達関数に基づいて、血流動態に関する情報を算出する血流情報算出手段と
を備える血流動態解析装置。
【請求項２】
前記第１の時刻は、前記画像データの収集開始時刻、または前記収集開始時刻と前記第２の時刻の間の所定時刻であり、前記第２の時刻は、前記動脈に対応する前記時間濃度曲線が上昇し始める時刻である請求項１に記載の血流動態解析装置。
【請求項３】
前記第１の集合は、１次元に配列されたＬ個の異なる時刻での前記組織における前記造影剤の濃度に対応する値の集合に対して、時間軸方向の前方にＨ個の第１の要素が追加された集合であり、前記第２の集合は、２次元に配列されたＬ個の異なる時刻での前記動脈における前記造影剤の濃度に対応する値の集合に対して、前記時間軸方向の後方にＨ個の第２の要素が追加された集合である請求項１又は２に記載の血流動態解析装置。
【請求項４】
前記第１の要素は、０及び前記組織における前記造影剤の濃度に対応する値の少なくとも一方であり、前記第２の要素は、０及び前記動脈における前記造影剤の濃度に対応する値の少なくとも一方である請求項３に記載の血流動態解析装置。
【請求項５】
前記血流情報算出手段は、前記伝達関数の最大値を検索し、検索した前記伝達関数の最大値からピーク範囲を算出し、前記伝達関数の最大値および前記ピーク範囲に基づいて、前記血流動態に関する情報を算出する請求項１乃至４のいずれか１項に記載の血流動態解析装置。
【請求項６】
コンピュータを、
医用画像撮影装置により得られた複数の時相における被検体の造影画像データに基づいて、前記複数の時相間における各画素の画素値の変化を計測して組織における造影剤の濃度に対応する量の時間変化を表す時間濃度曲線を算出する組織の時間濃度曲線算出手段、
前記組織に対応する前記時間濃度曲線から動脈領域を設定し、動脈に対応する時間濃度曲線を算出する動脈の時間濃度曲線算出手段、
前記組織に対応する前記時間濃度曲線を、第１の時刻を基準にして離散化し、異なる時刻での前記組織における前記造影剤の濃度に対応する値を含む複数の要素を１次元に配列した第１の集合を構成する第１の構成手段、
前記動脈に対応する前記時間濃度曲線を、前記第１の時刻より後方の第２の時刻を基準にして離散化し、異なる時刻での前記動脈における前記造影剤の濃度に対応する値を含む複数の要素を２次元に配列した第２の集合を構成する第２の構成手段、
前記第１の集合と前記第２の集合の逆畳み込み積分を行い、前記組織の伝達関数を算出する逆畳み込み積分手段、及び
前記伝達関数に基づいて、血流動態に関する情報を算出する血流情報算出手段、
として機能させる血流動態解析プログラム。
【請求項７】
被検体の複数のフレームの画像データに基づいて、組織に対応する画素において画素値の濃度に応じた量の前記複数のフレーム間における第１の変化を算出する一方、流体に対応する画素値の濃度に応じた量の前記複数のフレーム間における第２の変化を算出する濃度変化算出手段と、
前記第１の変化に基づいて前記第１の変化の第１の点に対応する要素を最初の要素とする解析対象となる数の第１の複数の要素を作成する一方、前記第２の変化に基づいて前記第２の変化の前記第１の点よりも後の第２の点に対応する要素を最初の要素とする前記解析対象となる数の第２の複数の要素を作成し、前記第１の複数の要素及び前記第２の複数の要素に基づいて逆畳み込み積分処理を伴って流体の動態に関する情報を算出する流体情報算出手段と、
を備える流体解析装置。
【請求項８】
前記濃度変化算出手段は、複数のフレームの造影画像データに基づいて、造影剤の濃度に応じた量の時間変化を前記第１及び第２の変化として算出するように構成される請求項７に記載の流体解析装置。
【請求項９】
前記濃度変化算出手段は、複数のフレームの非造影画像データに基づいて、前記第１及び第２の変化を算出するように構成される請求項７に記載の流体解析装置。
【請求項１０】
前記濃度変化算出手段は、磁気共鳴イメージング装置により共通のラベリングパルスからイメージングデータの収集タイミングまでの時間を変えてダイナミック収集された複数のフレームの非造影画像データに基づいて、前記第１及び第２の変化を算出するように構成される請求項９に記載の流体解析装置。
【請求項１１】
前記濃度変化算出手段は、磁気共鳴イメージング装置により反転回復パルス又は飽和パルスの印加とイメージングデータの収集とを繰り返し、かつ前記反転回復パルス又は飽和パルスの印加タイミングから前記イメージングデータの収集タイミングまでの時間を変えて収集された複数のフレームの非造影画像データに基づいて、前記第１及び第２の変化を算出するように構成される請求項９に記載の流体解析装置。
【請求項１２】
前記濃度変化算出手段は、磁気共鳴イメージング装置により反転回復パルス又は飽和パルスの印加とイメージングデータの収集とを繰り返し、かつ前記反転回復パルス又は飽和パルスの印加領域を変えて収集された複数のフレームの非造影画像データに基づいて、前記第１及び第２の変化を算出するように構成される請求項９に記載の流体解析装置。
【請求項１３】
コンピュータを、
被検体の複数のフレームの画像データに基づいて、組織に対応する画素において画素値の濃度に応じた量の前記複数のフレーム間における第１の変化を算出する一方、流体に対応する画素値の濃度に応じた量の前記複数のフレーム間における第２の変化を算出する濃度変化算出手段、及び
前記第１の変化に基づいて前記第１の変化の第１の点に対応する要素を最初の要素とする解析対象となる数の第１の複数の要素を作成する一方、前記第２の変化に基づいて前記第２の変化の前記第１の点よりも後の第２の点に対応する要素を最初の要素とする前記解析対象となる数の第２の複数の要素を作成し、前記第１の複数の要素及び前記第２の複数の要素に基づいて逆畳み込み積分処理を伴って流体の動態に関する情報を算出する流体情報算出手段、
として機能させる流体解析プログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【公開番号】特開２０１１−１３１０４１（Ｐ２０１１−１３１０４１Ａ）
【公開日】平成２３年７月７日（２０１１．７．７）
【国際特許分類】
【出願番号】特願２０１０−２００４３６（Ｐ２０１０−２００４３６）
【出願日】平成２２年９月８日（２０１０．９．８）
【出願人】（０００００３０７８）株式会社東芝 (54,554)
【出願人】（５９４１６４５４２）東芝メディカルシステムズ株式会社 (4,066)
【Ｆターム（参考）】
【公開日】平成２３年７月７日（２０１１．７．７）
【国際特許分類】
【出願日】平成２２年９月８日（２０１０．９．８）
【出願人】（０００００３０７８）株式会社東芝 (54,554)
【出願人】（５９４１６４５４２）東芝メディカルシステムズ株式会社 (4,066)
【Ｆターム（参考）】
[ Back to top ]