WO2004024001A1

WO2004024001A1 - 血流動態解析装置とその方法、画像診断装置

Info

Publication number: WO2004024001A1
Application number: PCT/JP2003/011700
Authority: WO
Inventors: Yasuo Omi; Osamu Miyazaki; Masahiko Aoyagi
Original assignee: Hitachi Medical Corporation
Priority date: 2002-09-12
Filing date: 2003-09-12
Publication date: 2004-03-25
Also published as: EP1537824A1; JP2004097665A; CN1681438A; US20060098852A1; US7436990B2; EP1537824A4; CN100377691C; JP4013128B2

Abstract

本発明は、ダイナミック撮影した断層像に基づいて流入動脈の時間−濃度曲線及び各組織における時間−濃度曲線をフーリエ変換し、流入動脈の時間−濃度曲線のフーリエ変換から逆フィルタを算出し、この逆フィルタを各組織における時間−濃度曲線のフーリエ変換に掛け合わせることで各組織の伝達関数を求め、このようにして求めた各組織の伝達関数を使用して生体機能情報を算出する。これにより、コンピュータ断層診断装置の提供する断層像から解析対象とする器官の生体機能情報を得る場合において、低い造影レートで定量性のよい生体機能情報を得ることができ、特にできる限り短い演算時間で生体機能情報を得ることが可能となる。

Description

明細書血流動態解析装置とその方法、画像診断装置技術分野

本発明は、血流動態解析装置とその方法及び画像診断装置に係り、特に X線や電磁波を用いて被検体の時間変化に依存する生 ί«能情報を得ることが可能なコンピュータ断層診断装置の提供する断層像から血流動態解析装置等の画像解析を行うことに対して有用な技術に関する。背景技術

従来の血流動態解析装置は、例えば、 X線 CT装置を用いてダイナミック撮影を行い、その撮影時に、ョード系の造影剤を患者に注入し、その注入された造影剤の濃度の時間変化の情報を得ている。その時間変化から時間一濃度曲線を得ることができ、この時間一濃度曲線の形状を解析することで、血流動態などの生体機能情報を得ることができる。その血流動態を解析する代表的なァルゴリズムは、ガンマ関数解析法、最大勾配法、デコンポリューション（deconvolution) 法がある。

まず、ガンマ関数解析法では、前記時間一濃度曲線をガンマ関数で近似し、その近似曲線のピーク値や曲線下面積から血流情報を算出する。次に、最大勾配法では、各組織における時間一濃度曲線の傾きの最大値を動脈入力関数における CT 値上昇の最大値で除算することで血流量を算出する。

ところが、ガンマ関数解析法及ぴ最大勾配法は、 1 秒当たりに注入する造影剤の量（造影レート）力 S 8〜10mlZ秒程度必要で、造影剤が注入される患者にとつて体力的な負担が大きいという問題がある。さらにガンマ関数解析法では定性評価はできるが、定量評価ができないと問題がある。

そこで、前記造影レートを低くでき、かつ定量評価が可能な血流解析方法が求められた結果、その一つとして、デコンボリューシヨン法が提案された（公知文献参照）このデコンポリューション法では、動脈入力関数と組織残留関数をデコンポリューション演算することでインパルス残留関数を求め、その求めたインパルス残留関数のピーク値や曲線下面積から血流情報を算出する。このデコンポリユーション法の利点は、 3〜5mlZ秒程度の低い造影レートで検査可能であるため、上記ガンマ関数解析法や最大勾配法と比較して半分程度の造影剤の注入速度でよい。

公知文献： L.Ostergaard etc.:High Resolution Measurement of Cerebral Blood Flow using Intravascular Tracer Bolus Passages: 1996； Magnetic Resonance in Medicine Vol.36:P.715-725)

しかしながら、デコンボリューシヨン法は、動脈入力関数と組織残留関数とからィンパルス残留関数へのデータ変換演算に際し複数回の積分をする必要があるため、その演算処理時間を要する。

本発明の目的は、短時間での血流動態解析を可能とすることにある。発明の開示

上記目的を達成するため、本発明は、断層像撮影装置により撮影された断層像を入力し、前記入力された断層像の各画素の時間変化の情報を表す時間一濃度曲線を求め、前記求められた断層像の各画素の時間—濃度曲線から流入動脈の時間 —濃度曲線を抽出し、前記抽出された流入動脈の時間一濃度曲線から逆フィルタ関数を求め、前記求められた逆フィルタ関数と前記断層像の各画素の時間一濃度曲線に基づき前記断層像の各画素の伝達関数を求め、前記求められた前記断層像の各画素の伝達関数を用いて血流動態解析像を求めることを特徴としている。これにより、例えば、流入動脈の時間一濃度曲線の逆フィルタと各組織における時間一濃度曲線との周波数空間上での乗算により各組織の伝達関数を求めることができ、従来のデコンボリューション法よりも演算時間の短縮ィ匕を図ることができる。

また、望ましい一実施例によれば、前記断層像内の最大連結画素領域を抽出し、前記抽出された最大違結画素領域を示すマスク画像を作成し、前記作成されたマスク画像を使用して前記入力される断層像のうちの当該マスク画像の領域以外のルームエアーや寝台、骨を含む不要な領域を除去することを特徴としている。これにより、血流解析に不要な領域を除去することで、解析対象となるデータ量が小さくなるので、上記生体機能情報算出手段等の演算時間の短縮ィ匕を図ることができる。

また、別の望ましい一実施例によれば、前記求められた断層像の各画素の時間一濃度曲線から流出静脈の時間一濃度曲線を求めるとともに、前記求められた流入動脈の時間一濃度曲線及び流出静脈の時間一濃度曲線のピーク値をそれぞれ求める手段と、前記流入動脈の時間一濃度曲線のピーク値が前記流出静脈の時間一濃度曲線のピーク値が一致するように流入動脈の時間一濃度曲線におけるパーシャルボリュームアベレージング（「PVA」と略記する）効果を補正することを特徴としている。 . 即ち、 X線 CT装置を用いて血流動態を検査する場合、 PVA効果により動脈を含むポクセルの CT値は本来の値よりも低下するため、解析結果の定量性が劣化するので、それを捕正することで、鮮明な血流動態解析ができる。図面の簡単な説明

図 1は本発明に係る血流動態角军析装置のハードウエア構成例を示すプロック図である。図 2は本発明に係る血流動態解析装置の断層像の入力から機能画像の表示までを示すフローチャートである。図 3は不要領域を削除する方法を示すフロ一チャートである。図 4は本実施の形態における PVA効果の捕正方法を示す図である。図 5は本実施の形態における PVA効果の他の補正方法を示す図である。図 6は本実施の形態における血管径とピーク値の算出手順を示すフローチヤ一トである。図 7は本実施の形態における逆フィルタ法を説明するための概念図である。図 8は本実施の形態における逆フィルタ法を説明するための概念図である。図 9は本実施の形態における逆フィルタ演算の処理手順を示すフローチヤ一トである。図 10は生体機能情報の算出方法の説明図である。図 11は高周波制御フィルタの算出方法 1の説明図である。図 12は高周波制御フィルタの算出方法 2の説明図である。図 13は高周波制御フィルタの算出方法 3の説明図である。発明を実施するための最良の形態以下、添付図面に従って本発明に係る血流動態解析装置の好ましい実施の形態について詳説する。

図 1は本発明に係る血流動態解析装置のハードウエア構成例を示すプロック図である。同図に示すように、この血流動態解析装置は、主としてコンピュータ断層像撮影装置で撮影した断層像を入力する断層像入力装置 1と、画像解析などの諸演算を行う演算装置 2とから構成されている。

断層像入力装置 1と演算装置 2は必ずしも別個の装置である必要はなく、一体ィ匕していてもよい。

演算装置 2は、断層像入力装置 1から断層像の取り込みを行うインターフェイス（I/F) 3 と、断層像や演算結果を一時的に格納するメモリ 4と、諸演算を行う中央処理装置（CPU) 5 と、演算結果や機能画像を記録するハードディスクなどの記録媒体 6と、マウスやキーボードなどの外部入力装置 7と、処理結果を表示する表示装置 8とから構成され、これらは共通バス 9によって接続されている。図 2は本発明に係る血流動態解析装置の断層像の入力から機能画像の表示までを示すフローチャートである。

まず、断層像入力装置 1から血流動態解析を行う被検体の断層像を選択し、選択された断層像を演算装置 2に入力する（ステップ 201) 。入力された断層像はメモリ 4に一時格納されるカゝ、または記録媒体 6に保存される。断層像入力装置 1 と演算装置 2とが一体ィ匕している場合には、外部入力装置 7により血流動態解析を行う被検体の断層像を選択し、選択された断層像を記録媒体 6から読み出し、メモリ 4に一時格納する。

続いて、 CPU5 により、断層像に写り込んでいるルームエアーや寝台、骨等、生体機能情報の解析には不要な領域を取り除く（ステップ 202) 。不要領域除去の方法については後述する。

次に、 CPU5 により、断層像から画素ごとの時間変化の情報を表す時間一濃度曲線を取り出す（ステップ 203) 。その後、 CPU5により、時間一濃度曲線から第一循環成分を取り出す（ステップ 204) 。第一循環成分を取り出すアルゴリズムは、ガンマ関数フィッティングや指数関数による外揷等、任意の公知のアルゴリズムでよい。続いて、 CPU5 により、動脈入力関数に用いる動脈（以下、流入動脈と記述する）と PVA効果の補正等の定量ィ匕に用いる静脈（以下、流出静脈と記述する）を決定する (ステップ 205) 。流入動脈と流出静脈は、操作者が断層像を見ながら手動で指定してもよく、時間一濃度曲線のピーク値及ぴピーク時間に基づいて解析対象である器官に対する流入動脈と流出静脈を自動的に決定してもよい。流入動脈と流出静脈を自動的に決定する場合には、まず、流入動脈と流出静脈の中心画素を算出する。流入動脈と流出静脈の中心画素を自動選択するには、各組織における時間一濃度曲線の最大値と最小値の差（以下、 ATDC と記述する）やピーク時間、ピーク値の特性の違いを利用すればよい。例えば、 ATDC がある閾値以上の画素の中から、ピーク値が所定の閾値以上で、かつ最もピーク時間の小さレ、画素を選ぶことにより、流入動脈の中心画素を自動選択することができ、また、 A TDC がある閾値以上の画素の中からピーク値が所定の閾値以上で、カ^ つ最もピーク時間の大きい画素を選ぶことにより、流出静脈の中心画素を自動選択することができる。次に、流入動脈と流出静脈の中心画素周辺の画素に対して、ピーク値が閾値以上の画素とそうでない画素を分離し、ピーク値が閾値以上の画素のうち、中心画素を含んだ連結画素を抽出することにより、流入動脈領域と流出静脈領域が決定される。尚、中心画素のみを流入動脈と流出静脈としてもよい。次に、 CPU5により、 PVA効果の補正を行う（ステップ 206) 。尚、 PVA効果の補正方法については後述する。続いて、 CPU5 により、画素毎に逆フィルタ演算を行い、生体機能情報を算出する (ステップ 207) 。尚、逆フィルタ演算については後述する。

画素毎に解析された生体機能情報をマッビングすることで機能画像が得られ (ステップ 208)、表示装置 8により機能画像が表示される（ステップ 209)。尚、ステップ 204は、ステップ 206とステップ 207の間で実行してもよい。

次に、不要領域を削除する方法について説明する。

図 3は不要領域を削除する方法を示すフローチャートである。同図に示すように、まず、 CPU5により、図 2のステップ 201で入力した断層像を 2値ィ匕する（ステツプ 301) 。 2値化処理は、例えば、断層像において、閾値以上の画素を 1に、閾値未満の画素を 0にそれぞれ置き換えることで達成される。なお、閾値は、解祈対象となる生体組織とルームエアーが分離できるような任意の値でよい。例えば、入力した断層像が CT像であるならば、閾値を一200程度に設定すればよレ、。続いて、 CPU5により、ステップ 301で 2値ィ匕された画像をラベリング処理する（ステップ 302) 。ラベリング処理は、つながつている全ての画素（連結成分）に同じ番号（ラベル）を付けて、異なる連結成分には別の番号を付けることで達成される。尚、隣接する画素が互いに 1の場合には、これらの画素はつながっていることになる。ラベルとして付与する値は、例えば 50からの続き番号とする。次に、 CPU5により、最大連結成分を探索する（ステップ 303) 。最大連結成分の探索は、ラベリングした画像全体を走査して、ラベル値ごとに画素数を数える処理を行い、最も画素数が多いラベル値を有する連結成分を選び出すことで達成される。その後、 CPU5により、最大連結成分を抽出する（ステップ 304)。最大連結成分の抽出は、ステップ 303で選び出された最も画素数が多、ラベル値を有する画素を残し、それ以外の画素を 0に置き換えることで達成される。これにより 2値化処理の段階では画素値が 1に置き換えられる可能性のある寝台や点滴チューブ等の不要領域の画素値が 0に置き換えられる。

次に、 CPU5 により、最大連結成分の輪郭を追跡し、最大連結成分の最外周の輪郭線を抽出する（ステップ 305) 。輪郭線の追跡は、ステップ 304で最大連結成分のみが残された画像の左上隅画素から横に走査し、最初に出会つた非ゼ口のラベル値を有する画素を開始点として反時計回り方向に輪郭を追跡し、開始点に戻ったら追跡を終了することで達成される。上記輪郭線上の画素にはラベル値とは異なる値、例えば 1を代入しておく。

続いて、 CPU5により、ステップ 305で求めた輪郭線の内部を 1で塗りつぶし、マスク画像を作成する（ステップ 306) 。輪郭線の内部を塗りつぶすには、従来から知られている、閉領域内部の 1点を開始点として閉領域内部を塗りつぶす処理であるシードフィルアルゴリズムを適用すればよい。本実施の形態では、閉領域はステップ 305で求めた輪郭線である。また、閉領域内部の点は、ステップ 304 で求めた最大連結成分のラベル値の付いた画素である。このラペル値の付いた画素を 1つ検出し、その点を開始点としてシードフィル処理を行えば、輪郭線内部を塗りつぶすことができる。次に、 CPU5により、入力画像（断層像）とステップ 306で作成したマスク画像とを乗算し、マスク領域外の画素を削除する（ステップ 307) 。具体的には、マスク画像を 1画素ずつ調べ、ゼロの画素を検出したら、その画素の座標に対応する原画像の画素を最低輝度に置き換える処理をする。これを全画素について行えば、マスク領域外の画素が入力画像から削除される。

次に、 PVA効果の補正方法について説明する。

図 4は本実施の形態における PVA効果の捕正方法を示す図である。例えば、脳の生体機能情報を解析する場合、一般に流入動脈には前大脳動脈や中大脳動脈を、流出静脈には上矢状静脈洞を用いる。 Lapin らの研究 (Journal of Computer Assisted Tomography 1993;Vol.l7： P.108— 114)によれば、 PVA効果により動脈を含むポクセルの信号値が低下するため、流入動脈である前大脳動脈や中大脳動脈の時間一濃度曲線は本来の値よりも小さくなる。また Lee らの研究 (American Journal of Neuroradiology 1999; Vol.20: P.63— 73)によれば、直径 1.73mm以上の血管は PVA効果の影響を受けないことから、流出静脈である上矢状静脈洞の時間 —濃度曲線は PVA効果の影響を受けていない。

よって流入動脈と流出静脈の時間濃度曲線は、一般に図 4 (a) のようになり、流入動脈の時間濃度曲線のピーク値 Paは、流出静脈のピーク値 Pvよりも PVA 効果による信号値の低下の分だけ小さな値になる。理想的には血管内の造影剤濃度は血管径の大小に依らず一定であるので、 PVA効果の影響が無ければ信号値のピーク値は血管径の大小によらず一定である。

従って、 .流入動脈のピーク値を Pa、流出静脈のピーク値を Pvとすると、図 4 (b) に示すように流入動脈の時間一濃度曲線を PvZPa倍すれば、流入動脈における PVA効果を捕正できる。 PVA効果を捕正する前の時間 tにおける流入動脈の時間一濃度曲線を AIFp_Va(t) 、 PVA効果を捕正した彼の時間 tにおける流入動脈の時間一濃度曲線を AIF(t)とすると、両者の間に次式が成り立つように補正すればよい。

AIF{t) = AIFpva(t)x— . . . (i)

Pa 図 5は本実施の形態における PVA効果の他の捕正方法を示す図である。過去の臨床画像の中に、これから実施する検査と全く同一の撮影条件、造影条件で撮影した例があるならば、この過去の臨床画像を用いて PVA効果を捕正することも可能である。

例えば、図 5 (a) に示す臨床画像において、断層像内に含まれる複数の血管 1 〜5に対し、血管径 R1〜！5とその血管の時間一濃度曲線のピーク値を算出する。これを図表に表すと、図 5 (b) となり、プロットしてグラフに表すと、図 5 (c) のようなカーブが得られる。尚、血管径とピーク値の算出方法については後述する。血管径が細いほど PVA効果の影響を強く受け、血管径がある値以上であれば PVA効果の影響を受けないため、図 5 (c) のカーブはある径 Rthまでは徐々に上昇し、 Rth以上の Rに対してはピーク値が平衡になる。尚、この平衡状態におけるピーク値を便宜上 Pt とし、流入動脈における血管径を Ra、流入動脈における時間一濃度曲線のピーク値を Paとする。

Raが Rthより小さい場合には、流入動脈の時間一濃度曲線を PthZPa倍すれば、流入動脈における PVA効果を補正できる。 PVA効果を捕正する前の時間 tにおける流入動脈の時間一濃度曲線を AEFpva(t)、 PVA効果を補正した後の時間 tにおける流入動脈の時間一濃度曲線を AIF(t)とすると、両者の間に次式が成り立つように補正すればよい。

AIF{t) = AIFpva{t)x— . . . (2)

Pa 次に、血管径とピーク値の算出方法について説明する。

図 6は本実施の形態における血管径とピーク値の算出手順を示すフローチヤ一トである。まず、 CPU5 により、各時間一濃度曲線のピーク時間とピーク値を算出する（ステップ 601) 。次に、 CPU5により、血管の中心画素を算出する（ステップ 602) 。中心画素を算出するには、ピーク値の特性を利用すればよく、ピーク値が最大の画素を選ぶことにより、血管の中心画素を自動選択することができる。続いて、 CPU5 により、血管の中心画素周辺の画素に対して、ピーク値が閾値以上の画素とそうでない画素を分離する（ステップ 603) 。その後、ステップ 603 で求めたピーク値が閾値以上の画素のうち、中心画素を含んだ連結画素を抽出することにより、血管領域が決定される（ステップ 604) 。次に、 CPU5により、血管領域内の画素数を算出し（ステップ 605) 、その算出した画素数から血管領域の面積を算出する（ステップ 606) 。一画素の面積は装置毎に一定値であるので、一画素あたりの面積にステップ 605で求められた画素数を掛けることで血管領域の面積 Sが得られる。

続いて、 CPU5により、血管径 Rを算出する（ステップ 60カ。断層像上での血管の形状は円形で近似して差し支えないので、血管径 Rと面積 Sには、次式が成り立ち、その式に従って血管径を算出すればよい。

、

R = (3) 次に、 CPU5により血管の時間一濃度曲線を算出する（ステップ 608)。ステツプ 203において既に各画素における時間一濃度曲線が算出されているので、血管領域の全ての画素における時間一濃度曲線の平均値を算出することで、血管の時間一濃度曲線が得られる。また、平均値をとるかわりに、ステップ 602で算出されている血管中心画素からの距離に応じて、任意に重み付けをしてもよい。

次に、 CPU5 により、血管の時間一濃度曲線のピーク値を算出する (ステップ 609) 。尚、ステップ 605からステップ 607までの処理と、ステップ 608からステップ 609までの処理は、処理の順序が逆になつてもよい。本実施例では、ピーク値の特性を利用して血管の抽出を行ったが、血管の抽出法はこれに限定されるものではなく、血管が正確に抽出できる方法であるならば任意の方法でよい。また、断層像を見ながら手動で血管領域を指定することにより、血管を抽出してもよい。次に、逆フィルタ法の概念について説明する。

図 7は本実施の形態における逆フィルタ法を説明するための概念図である。伝達関数 f の対象組織において、入力関数を On (り、出力関数を Cout(t) とおくと Cout(t) は Cin(t)と fを畳み込み積分したものになり、次式のように表される。

Cout t) = dn{t)- f(t - f)df (4) 造影剤を動脈から体内に注入する場合、 On(t)は流入動脈における時間一濃度曲線（動脈入力関数） AIF(t)に、 Cout(t) は組織における時間一濃度曲線 Cr(t) に置き換えられ、次式のように表される。 Cr{t) = )AIF{t)- f{t - f)df · · · (5)

0

(5) 式を周波数空間に変換すると、次式のようになる。

ί

F{Cr(t)}= F{AIF(t)}- F{f{t)} . . . (6) ここで、 F {AEF(t)} は AIF(t)のフーリエ変換を、 F {Cr(t) } は Cr(t) のフーリェ変換を、 F {f(t)} は f(t)のフーリエ変換をそれぞれ表す。

(5) 式と（6) 式より、流入動脈における時間一濃度曲線と組織における時間 —濃度曲線を逆畳み込み演算すれば伝達関数 f を求めることができるが、逆畳み込み演算を行う代わりに、流入動脈における時間一濃度曲線の逆フィルタを組織の時間一濃度曲線に掛けても等価な演算結果が得られる。

AIF(t)の逆フィルタを AIF— )とおくと、（6) 式より、 f(t)は次式のように表される。 f{t) =

· ■ · (フ）ここで F— ¹は、逆フーリエ変換を表す。

図 8は本実施の形態における逆フィルタ法を説明するための概念図である。血流動態解析を行う場合、被検者への侵襲性を抑えるため、たいていの場合、造影剤は静脈から体内へ注入する。この場合、造影剤は体内を循環した後、解析対象組織に流入する。入力関数を Cin (り、循環経路における伝達関数を、対象組織における伝達関数を ¾、対象組織からの出力関数を Cout(t) とおくと、 Cout(t) は、 Cin(t)と &を畳み込み積分したものを、更に f₂と畳み込み積分したものになるので、次式のようになる。

Cout{t) = f{|Cm(t) - - t^dt,)- f₂{t - t₂')dt₂' · · . (8)

0

Cin(t)との畳み込み積分は流入動脈における時間一濃度曲線（動脈入力関数） AIF(t)に、 Cout(t) は組織における時間—濃度曲線 Cr(t) に置き換えられ、次式のように表される。 Cr(t) = f₂{t - t₂')dt₂' · · · (9)

(9) 式は、前述した（5) 式と等価である。よって伝達関数 f₂は、次式のように表される。 2(t) = F-¹[E{A/E-¹(t)}- E{Cr(t)}] · · · (10) 次に、逆フィルタ演算について説明する。

図 9は本実施の形態における逆フィルタ演算の処理手順を示すフローチヤ一トである。同図に示すように、 CPU5 により、流入動脈の時間一濃度曲線をフーリェ変換する（ステップ 901) 。次に、 CPU5により、流入動脈の時間—濃度曲線のフーリエ変換から逆フィルタを算出する（ステップ 902) 。流入動脈の時間一濃度曲線のフーリエ変換の実部を F {AIF_R} 、虚部を F {AIF とすると、逆フィルタは、次式で表される。

F{AIF_R} iF{AIF_r}

[FIAIF +IAIF ] [F{AIF_R +{AIF_I}²^

次に、逆フィルタの高周波数成分を制御する（ステップ 903) 。ウィナーフィルタゃバターワースフィルタなどの高周波制御フィルタによって逆フィルタをフィルタリングすることにより、逆フィルタの高周波数成分が制御される。続いて、 CPU5により、各組織の時間一濃度曲線をフーリェ変換する（ステップ 904)。このフーリエ変換に対して、ステップ 902で求めた逆フィルタでフィルタリングする（ステップ 905) 。このフィルタリングされた糸且織の時間一濃度曲線のフーリエ変換を逆フーリエ変換する伝達関数が求められる（ステップ 906)。

次に、 CPU5により、ステップ 906で求められた伝達関数から、血流動態などの生体機能情報を算出する（ステップ 907) 。生体機能情報の算出方法については後述する。尚、ステップ 904からステップ 907までは、解析する組織の断層像上での画素数分だけ繰り返す。また、ステップ 901からステップ 903までと、ステツプ 904は、処理の順序は逆になつてもよい。

次に、伝達関数から生 ί«能情報を算出する方法について説明する。

図 10は生 ί«能情報の算出方法の説明図である。伝達関数を f(t)とすると、 f(t) は、例えば図 10に示すような形状になり、血流量（blood flow:「BF」と略記する）は、次式のように f(t)の最大値から求められる。

BF = f_max · · . (12)

また、血液量（blood volume:以下、 BVと記述する）は、次式のように f(t)の曲線下面積から求められる。

平均通過時間（meantransit time: 以下、 MTT と記述する）は、次式のように f(t) の幅から求められる。

ここで、 f maxは、 f(t)の最大値を表す。

続いて、高周波制御フィルタの決定方法については説明する。式（9) から伝達関数を算出するには、動脈における時間—濃度曲線と組織における時間一濃度曲線とを逆たたみこみ演算する必要がある。逆たたみこみ演算は大きく分けて、特異値分解法のように行列式を用いる方法と逆フィルタ法のようにフーリエ変換を用いる方法とが存在し、フーリェ変换を用レヽる方法は行例式を用ヽる方法に比べて演算時間を短縮できる。逆フィルタ法では動脈における時間一濃度曲線のフーリェ変換から算出した逆フイノレタで、組織における時間一濃度曲線のフーリェ変換を周波数空間上でフィルタリングするが、この演算は高周波数成分を強調する処理である。一般に生体学的信号成分は低周波数領域に、ノイズ成分は高周波数成分に現われるため、単純に逆フィルタ演算を行っただけではノイズを無用に強調してしまうため生 ί«能情報を精度よく得ることは難しい。このため逆フィルタに対して高周波数成分を制御する必要がある。

対象組織への動脈入力関数、すなわち流入動脈における時間一濃度曲線が図

11a のように実測されたとすると、この時間一濃度曲線の平均通過時間 (MTT_Arteiy) は、次式で求められる。 MTT = · · · (15)

ここで tは時間、 nは計測点数、 Ca(t)は時間—濃度曲線を表す。流入動脈における時間一濃度曲線のピークタイム（図 11中の PT) と、（16) 式から求められる平均通過時間から、図 11 (b) のような、ピークタイムを中心に平均通過時間の幅を持つ理想的な矩开多の入力関数 ideal_Ca(t)が求められる。この ideal_Ca(t)のフリエ変換を F{ideal— Ca(t)}とすると、 F{ideal_Ca(t)}は図 11 (c) のようになる。ここで有用な生態学的情報は低周波数成分として現われ解析に無意味なノィズ成分は低周波数成分に現われるので、 F{ideal— Ca(t)}に対して、例えば図 11 (c) に示すように角周波数 0から 2 Zkまでの成分のみをそのまま残し、 2π力、ら 4π までの成分を減弱し、 4πよりも高い周波数成分はカツトするようにフィルタの形状を決定する。この高周波制御フィルタは、次式で与えられる。

1 0 < 03 < ω ₅ filter、(o) = . . . (if.)

^o

0 ω <ω ここで、 ωは角周波数、 iosは制御開始周波数、 co eは制御終了周波数であり、 ideal—Ca(t)の幅、すなわち Ca(t)の平均通過時間を kとすれば、

coe =4 % /kである。よつて流入動脈における時間一濃度曲線から、高周波領域制御フィルタの形状を決定することができる。本実施の形態では、便宜上、 ω8=2π Zk coe =4π/]αとして説明したが、 cosと c eはこれに限られるものではなく、

(m≤n) として任意に設定してもよい。なお本実施の形態では、流入動脈における時間一濃度曲線から sと coeを算出したが、 cosと coeの算出方法はこれに限られるものではなく、図 12のように最大周波数から任意に設定してもよい。

高周波制御フィルタは、次式のようなフィルタを用いてもよレ、。 ' filter {ω) = filter 2{ω) + 2(ω)— / 1(ω)}

1 0 < ω < «2^ filterlifo) — co²j} ω2₃≤ω≤ 2_e

ω2 ≤ω

(17) 式のフィルタは、低周波数成分は強調しつつ高周波数成分は抑制するようなフィルタであり、 kは低周波数成分の強調度を表すパラメータである。ここで例えば w ls=0、 ω le= π o 2s=0、 ω 2^=3 π / k=0.75 と設定すれば図 13に示すような高周波制御フィルタが得られる。なお本実施の形態では、流入動脈における時間一濃度曲線から co ls、 ω 1β co 2s、 2e、を算出したが、算出方法はこれに限られるものではなく、最大周波数から任意に設定してもよい。また、上記高周波制御フィルタは、コンピュータ断層診断装置に備えられている画像表示部や計測の進行を表示部へグラフを表示してもよい。そのグラフには、マウスやトラックボールなどのボインティングデバイスを用いた周知のグラフィカル ·ユーザ 'インターフェースを用い、例えば表示されたグラフの形状を変形するようにして、上記高周波制御フィルタに与えられる各種パラメータの少なくとも 1つを変更できるようにしてもよい。各種パラメータは、高周波制御フィタの制御周波数、低域周波数領域の強調の度合いなどである。

これにより、各種パラメータを変更することで、例えば、高周波制御フィルタの制御周波数の帯域を拡げたり、低域周波数領域の強調の度合いを増加すれば、本来の信号を適宜に強調すれば、結果として信号成分が増加するから雑音成分が抑制されるための最適な高周波制御フィルタを提供することができる。よって、高精細な血流動態などの生体機能情報解析を行った画像を提供することができる。以上説明したように上記実施例によれば、コンピュータ断層診断装置によりダィナミック撮影された断層像から生能情報解析を行う場合、ガンマ関数解析？去ゃ最大勾配法よりも低い造影レートで、デコンポリューション法よりも短い演算時間で、生 ί«能情報^ fを行えるという効果がある。また、 PVA効果による定量性の悪ィヒを防止できるという効果もある。更に、操作者への負担を増加することなく不要領域の影響による無用な演算時間の増加や解析精度を防止することができるという効果もある。

Claims

6

求の

1. 断層像撮影手段により撮影された断層像を入力する入力手段と、

この入力手段により入力された断層像の各画素の時間変化の情報を表す時間一濃度曲線を求める第 1の演算手段と、

この第 1の演算手段により求められた断層像の各画素の時間一濃度曲線力ら流入動脈の時間—濃度曲線を抽出する第 2の演算手段と、

この第 2の演算手段により抽出された流入動脈の時間一濃度曲線から逆フィルタ関数を求める第 3の演算手段と、

この第 3の演算手段により求められた逆フィルタ関数と前記第 1の演算手段により求められた断層像の各画素の時間一濃度曲線に基づき前記断層像の各画素の伝達関数を求める第 4の演算手段と、

この第 4の演算手段により求められた前記断層像の各画素の伝達関数を用いて血流動態解析像を求める第 5の演算手段と、

を備えたことを特徴とする血流動態解析装置。

2. 前記第 3の演算手段は、前記求められた逆フィルタ関数の高周波成分を制御する高周波制御フィルタ手段を備えたことを特徴とする請求項 1に記載の血流動態解析装置。

3. 前記高周波制御フィルタ手段は、前記第 2の演算手段に抽出された流入動脈の時間一濃度曲線からフィルタ関数を求めるフィルタ関数演算手段を備えたことを特徴とする請求項 2に記載の血流動態解析装置。

4. 前記フィルタ関数演算手段は、前記第 2の演算手段に抽出された流入動脈の時間一濃度曲線をフーリエ変換し、そのフーリエ変換された流入動脈の時間一濃度曲線から前記高周波制御フィルタ手段の制御開始周波数と制御終了周波数を設定し、それら設定された制御開始周波数と制御終了周波数に基づき前記高周波制御フィルタ手段のフィルタ関数を求めることを特徴とする請求項 3に記載の血流動態解析装置。

5. 前記制御開始周波数と制御終了周波数は、前記フーリエ変換された流入動脈の時間一濃度曲線の最大周波数に基づいて求めることを特徴とする請求項 4に記載の血流動態解析装置。

前記高周波制御フィルタ手段は、その高周波制御フィルタの低周波数部分を強調する度合ヽと前記高周波制御フィルタの帯域との少なくとも一方のパラメータを設定するパラメータ設定手段を備えたことを特徴とする請求項 2に記載の血流動態解析装置。

前記パラメータ設定手段は、前記高周波制御フィルタ手段のフィルタ関数を表示する表示手段と、この表示手段により表示されたフィルタ関数を変更する手段と、を備えたことを特徴とする請求項 6に記載の血流動態解前記変更手段は、前記表示手段に表示された高周波制御フィルタ関数の形状をグラフィック ·ユーザ 'インターフェースにより変更することを特徴とする請求項 7に記載の血流動態解析装置。

前記入力手段により入力された断層像の内の最大連結画素領域を抽出する手段と、この抽出手段により抽出された最大連結画素領域に基づき前記入力手段から入力される断層像のうち血流動態解析に不要な領域を除去する手段と、を備えたことを特徴とする請求項 1に記載の血流動態解析

10. 前記不要領域はルームエアーや寝台、骨を含むことを特徴とする請求項 9に記載の血流動態解析装置。

11. 前記入力手段により入力された断層像の各画素の時間一濃度曲線から流出静脈の時間一濃度曲線を求めるとともに、その求められた流入動脈の時間一濃度曲線及び流出静脈の時間一濃度曲線のピーク値をそれぞれ求める手段と、前記求められた流入動脈の時間一濃度曲線のピーク値と前記流出静脈の時間一濃度曲線のピーク値とがー致するように流入動脈の時間一濃度曲線におけるパーシャルボリュームアベレージング効果を補正する補正手段と、を備えたことを特徴とする請求項 1に記載の血流動態解

12. 断層像撮影装置により撮影された断層像を入力する入力」

、力された断層像の各画素の時間変化の情報を表す時間一濃度曲線を求める第 1の演算二

この第 1の演算ステップにより求められた断層像の各画素の時間一濃度曲線から流入動脈の時間一濃度曲線を抽出する第 2の演算ステップと、この第 2の演算ステツプにより抽出された流入動脈の時間—濃度曲線から逆フィルタ関数を求める第 3の演算ステップと、

この第 3の演算ステップにより求められた逆フィルタ関数と前記第 1の演算手段により求められた断層像の各画素の時間一濃度曲線に基づき前記断層像の各画素の伝達関数を求める第 4の演算ステップと、

この第 4の演算ステツプにより求められた前記断層像の各画素の伝達関数を用いて血流動態解析像を求める第 5の演算ステップと、

を備えたことを特徴とする血流動態解析方法。

13. 前記第 3の演算ステップは、前記求められた逆フィルタ関数の高周波成分を制御する高周波制御フィルタステップを備えたことを特徴とする請求項 12に記載の血流動態解析方法。

14. 前記高周波制御フィルタステップは、前記第 2の演算ステップに抽出された流入動脈の時間一濃度曲線からフィルタ関数を求めるフィルタ関数演算ステップを備えたことを特徴とする請求項 13 に記載の血流動態解析方法。

15. 前記フィルタ関数演算ステップは、前記第 2の演算ステップに抽出された流入動脈'の時間一濃度曲線をフーリエ変換し、そのフーリエ変換された流入動脈の時間一濃度曲線から前記高周波制御フィルタ手段の制御開始周波数と制御終了周波数を設定し、それら設定された制御開始周波数と制御終了周波数に基づき前記高周波制御フィルタ手段のフィルタ関数を求めることを特徴とする請求項 14に記載の血流動態解析方法。

16. 前記制御開始周波数と制御終了周波数は、前記フーリエ変換された流入動脈の時間一濃度曲線の最大周波数に基づいて求めることを特徴とする請求項 15に記載の血流動態解析方法。

17. 前記高周波制御フィルタステップは、その高周波制御フィルタの低周波数部分を強調する度合いと前記高周波制御フィルタの帯域との少なくと 9

も一方のパラメータを設定するパラメータ設定ステップを備えたことを特徴とする請求項 13に記載の血流動態解析方法。

18. 前記パラメータ設定ステップは、前記高周波制御フィルタ手段のフィルタ関数を表示器に表示する表示ステップと、この表示ステツプにより表示されたフィルタ関数を変更するステップと、を備えたことを特徴とする請求項 17に記載の血流動態解析方法。

19. 前記変更ステップは、前記表示ステップに表示された高周波制御フィルタ関数の形状をグラフィック ·ユーザ 'インターフェースにより変更することを特徴とする請求項 18に記載の血流動態解析方法。

20. 前記入力ステップにより入力された断層像の内の最大連結画素領域を抽出するステップと、この抽出ステップにより抽出された最大連結画素領域に基づき前記入力ステップから入力される断層像のうち血流動態解析に不要な領域を除去するステップと、を備えたことを特徴とする請求項 12 に記載の血流動態解析方法。

21. 前記不要領域はルームエアーや寝台、骨を含むことを特徴とする請求項 20に記載の血流動態解析方法。

22. 前記入力手段により入力された断層像の各画素の時間一濃度曲線から流出静脈の時間一濃度曲線を求めるとともに、その求められた流入動脈の時間一濃度曲線及び流出静脈の時間一濃度曲線のピーク値をそれぞれ求めるステップと、前記求められた流入動脈の時間一濃度曲線のピーク値と前記流出静脈の時間—濃度曲線のピーク値とがー致するように流入動脈の時間一濃度曲線におけるパーシャルボリユームァベレージング効果を捕正する補正ステップと、を備えたことを特徴とする請求項' 12に記載の血流動態解析方法。

23. 生体から断層像を撮影可能な画像診断装置において、

前記撮影された断層像を入力する入力手段と、

この入力手段より入力された断層像の各画素の時間変化の情報を表す時間一濃度曲線を求める第 1の演算手段と、

この第 1の演算手段により求められた断層像の各画素の時間一濃度曲線から流入動脈の時間一濃度曲線を抽出する第 2の演算手段と、

この第 3の演算手段により求められた逆フィルタ関数と前記第 1の演算手段により求められた断層像の各画素の時間—濃度曲線に基づき前記断層像の各画素の伝達関数を求める第 4の演算手段と、

を備えたことを特徴とする画像診断装置。

24. 前記第 3の演算手段は、前記求められた逆フィルタ関数の高周波成分を制御する高周波制御フィルタ手段を備えたことを特徴とする請求項 23 に

25. 前記高周波制御フィルタ手段は、前記第 2の演算手段に抽出された流入動脈の時間一濃度曲線からフィルタ関数を求めるフィルタ関数演算手段を備えたことを特徴とする請求項 24に記載の画像診断装置。

26. 前記フィルタ関数演算手段は、前記第 2の演算手段に抽出された流入動脈の時間一濃度曲線をフーリエ変換し、そのフーリエ変換された流入動脈の時間一濃度曲線から前記高周波制御フィルタ手段の制御開始周波数と制御終了周波数を設定し、それら設定された制御開始周波数と制御終了周波数に基づき前記高周波制御フィルタ手段のフィルタ関数を求めることを特徴とする請求項 25に記載の画像診断装置。

27. 前記制御開始周波数と制御終了周波数は、前記フーリエ変換された流入動脈の時間一濃度曲線の最大周波数に基づいて求めることを特徴とする請求項 26に記載の画像診断装置。

28. 前記高周波制御フィルタ手段は、その高周波制御フィルタの低周波数部分を強調する度合いと前記高周波制御フィルタの帯域との少なくとも一方のパラメータを設定するパラメータ設定手段を備えたことを特徴とする請求項 24に記載の画像診断装置。

29. 前記パラメータ設定手段は、前記高周波制御フィルタ手段のフィルタ関 2

数を表示する表示手段と、この表示手段により表示されたフィルタ関数を変更する手段と、を備えたことを特徴とする請求項 28 に記載の画像診断

30. 前記変更手段は、前記表示手段に表示された高周波制御フィルタ関数の形状をグラフィック ·ユーザ 'インターフェースにより変更することを特徵とする請求項 29に曾己載の画像診断装置。