(57)【特許請求の範囲】
【請求項１】
断層像撮影手段により撮影された断層像を入力する入力手段と、
この入力手段により入力された断層像の各画素の時間変化の情報を表す時間−濃度曲線を求める第１の演算手段と、
この第１の演算手段により求められた断層像の各画素の時間−濃度曲線から流入動脈の時間−濃度曲線を抽出する第２の演算手段と、
この第２の演算手段により抽出された流入動脈の時間−濃度曲線から逆フィルタ関数を求める第３の演算手段と、
この第３の演算手段により求められた逆フィルタ関数と前記第１の演算手段により求められた断層像の各画素の時間−濃度曲線に基づき前記断層像の各画素の伝達関数を求める第４の演算手段と、
この第４の演算手段により求められた前記断層像の各画素の伝達関数を用いて血流動態解析像を求める第５の演算手段と、
を備えたことを特徴とする血流動態解析装置。
【請求項２】
前記第３の演算手段は、前記求められた逆フィルタ関数の高周波成分を制御する高周波制御フィルタ手段を備えたことを特徴とする請求項１に記載の血流動態解析装置。
【請求項３】
前記高周波制御フィルタ手段は、前記第２の演算手段に抽出された流入動脈の時間―濃度曲線からフィルタ関数を求めるフィルタ関数演算手段を備えたことを特徴とする請求項２に記載の血流動態解析装置。
【請求項４】
前記入力手段により入力された断層像のうちの最大連結画素領域を抽出する手段と、この抽出手段により抽出された最大連結画素領域に基づき前記入力手段から入力される断層像のうち血流動態解析に不要な領域を除去する手段と、を備えたことを特徴とする請求項１に記載の血流動態解析装置。
【請求項５】
前記不要領域は、ルームエアー、寝台、及び骨を含むことを特徴とする請求項４に記載の血流動態解析装置。
【請求項６】
前記入力手段により入力された断層像の各画素の時間−濃度曲線から流出静脈の時間−濃度曲線を求めるとともに、その求められた流入動脈の時間−濃度曲線及び前記流出静脈の時間−濃度曲線のピーク値をそれぞれ求める手段と、前記求められた流入動脈の時間−濃度曲線のピーク値と前記流出静脈の時間―濃度曲線のピーク値とが一致するように流入動脈の時間−濃度曲線におけるパーシャルボリュームアベレージング効果を補正する補正手段と、を備えたことを特徴とする請求項１に記載の血流動態解析装置。
【請求項７】
断層像撮影装置により撮影された断層像を入力する入力ステップと、
この入力ステップにより入力された断層像の各画素の時間変化の情報を表す時間−濃度曲線を求める第１の演算ステップと、
この第１の演算ステップにより求められた断層像の各画素の時間−濃度曲線から流入動脈の時間−濃度曲線を抽出する第２の演算ステップと、
この第２の演算ステップにより抽出された流入動脈の時間−濃度曲線から逆フィルタ関数を求める第３の演算ステップと、
この第３の演算ステップにより求められた逆フィルタ関数と前記第１の演算ステップにより求められた断層像の各画素の時間−濃度曲線に基づき前記断層像の各画素の伝達関数を求める第４の演算ステップと、
この第４の演算ステップにより求められた前記断層像の各画素の伝達関数を用いて血流動態解析像を求める第５の演算ステップと、
を含むことを特徴とする血流動態解析方法。
【請求項８】
前記第３の演算ステップは、前記求められた逆フィルタ関数の高周波成分を制御する高周波制御フィルタステップを含むことを特徴とする請求項７に記載の血流動態解析方法。
【請求項９】
前記高周波制御フィルタステップは、前記第２の演算ステップに抽出された流入動脈の時間―濃度曲線からフィルタ関数を求めるフィルタ関数演算ステップを含むことを特徴とする請求項８に記載の血流動態解析方法。
【請求項１０】
前記入力ステップにより入力された断層像のうちの最大連結画素領域を抽出するステップと、この抽出ステップにより抽出された最大連結画素領域に基づき前記入力ステップから入力される断層像のうち血流動態解析に不要な領域を除去するステップと、を含むことを特徴とする請求項７に記載の血流動態解析方法。
【請求項１１】
前記不要領域は、ルームエアー、寝台、及び骨を含むことを特徴とする請求項１０に記載の血流動態解析方法。
【請求項１２】
前記入力ステップにより入力された断層像の各画素の時間−濃度曲線から流出静脈の時間−濃度曲線を求めるとともに、その求められた流入動脈の時間−濃度曲線及び前記流出静脈の時間−濃度曲線のピーク値をそれぞれ求めるステップと、前記求められた流入動脈の時間−濃度曲線のピーク値と前記流出静脈の時間―濃度曲線のピーク値とが一致するように流入動脈の時間−濃度曲線におけるパーシャルボリュームアベレージング効果を補正する補正ステップと、を含むことを特徴とする請求項７に記載の血流動態解析方法。
【請求項１３】
生体から断層像を撮影可能な画像診断装置において、
前記撮影された断層像を入力する入力手段と、
この入力手段により入力された断層像の各画素の時間変化の情報を表す時間−濃度曲線を求める第１の演算手段と、
この第１の演算手段により求められた断層像の各画素の時間−濃度曲線から流入動脈の時間−濃度曲線を抽出する第２の演算手段と、
この第２の演算手段により抽出された流入動脈の時間−濃度曲線から逆フィルタ関数を求める第３の演算手段と、
この第３の演算手段により求められた逆フィルタ関数と前記第１の演算手段により求められた断層像の各画素の時間−濃度曲線に基づき前記断層像の各画素の伝達関数を求める第４の演算手段と、
この第４の演算手段により求められた前記断層像の各画素の伝達関数を用いて血流動態解析像を求める第５の演算手段と、
を備えたことを特徴とする画像診断装置。
【請求項１４】
前記第３の演算手段は、前記求められた逆フィルタ関数の高周波成分を制御する高周波制御フィルタ手段を備えたことを特徴とする請求項１３に記載の画像診断装置。
【請求項１５】
前記高周波制御フィルタ手段は、前記第２の演算手段に抽出された流入動脈の時間―濃度曲線からフィルタ関数を求めるフィルタ関数演算手段を備えたことを特徴とする請求項１４に記載の画像診断装置。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は血流動態解析装置に係り、特にＸ線や電磁波を用いて被検体の時間変化に依存する生体機能情報を得ることが可能なコンピュータ断層診断装置の提供する断層像から生体機能情報解析等の画像解析を行うことに対して有用な技術に関する。
【０００２】
【従来の技術】
従来、ダイナミック撮影を行うコンピュータ断層診断装置として、陽電子放射断層撮影装置（Positron Emission Tomography:PET）や単光子放射ＣＴ（Single Photon Emission CT:SPECT)などの核医学装置がある。核医学装置による血流動態計測は、放射性核種を被検体内に注入し、核種からの放射線をシンチレーションカメラで計測すると、横断面の放射性核種分布像を得られ、この放射性核種分布像を解析することで、器官の生体機能情報を解析していた。
【０００３】
他のダイナミック撮影を行うコンピュータ断層診断装置として、Ｘ線ＣＴ装置がある。単純Ｘ線ＣＴ像では病変の超早期相の診断が困難である場合には、Ｘ線ＣＴ装置によるダイナミック撮影を行い、この撮影時にコントラストを付けるために主にヨード系の造影剤を用いていた。
【０００４】
造影剤を注入してダイナミックスキャンを行うと、造影剤濃度とＣＴ値は比例関係にあることから、撮影断面の時間変化の情報が得られる。撮影断面の時間変化から各組織のＣＴ値の時間変化である時間−濃度曲線が得られ、この時間−濃度曲線を解析することで器官の生体機能情報を解析していた（例えば非特許文献１参照。）。
【０００５】
更に他のダイナミック撮影を行うコンピュータ断層診断装置として、ＭＲ装置がある。高磁場ＭＲ装置では非造影の診断方法もあるが、高磁場ＭＲ装置でコントラストの良い画像を得たい場合や低磁場ＭＲ装置においては、主にガドリニウム系の造影剤を用いていた。ＭＲ装置においても撮影断面の時間変化から各組織のＭＲ値の相対的信号強度変化率である時間−濃度曲線が得られる。この時間−濃度曲線を解析することで器官の時間変化に依存する生体機能情報を解析していた（例えば非特許文献２参照。）。
【０００６】
一方、血流動態を解析するアルゴリズムとして、ガンマ関数解析法がある。ガンマ関数解析法では、時間−濃度曲線をガンマ関数でフィッティングし、フィッティングカーブのピーク値や曲線下面積から血流情報を算出していた（例えば非特許文献３参照。）。
【０００７】
他の血流動態を解析するアルゴリズムとして、最大勾配法（maximum slope 法）がある。最大勾配法では、各組織における時間−濃度曲線の傾きの最大値を動脈入力関数におけるＣＴ値上昇の最大値で割ることで血流量を算出していた（例えば非特許文献４参照。）。
【０００８】
更に他の血流動態を解析するアルゴリズムとして、deconvolution(デコンボリューション）法がある。deconvolution 法では、動脈入力関数と組織残留関数をdeconvolution 演算することでインパルス残留関数を求め、インパルス残留関数のピーク値や曲線下面積から血流情報を算出していた（例えば非特許文献５参照。）。
【０００９】
【非特許文献１】
J.P.Muizelaar etc.:A New Method for Quantitative Regional
Cerebral Blood
Volume Measurements Using Computed Tomography.:Stroke
1997;28:1998-2005
【００１０】
【非特許文献２】
井川房雄  他：拡散強調画像と潅流画像による急性期脳梗塞の診断の意義：日磁医誌  第２０巻５号２０３−２０９頁（２０００年）
【００１１】
【非特許文献３】
L.Axel:Cerebral Blood Flow Determination by Rapid-Sequence
Computed Tomography:Radiology 1980;137:679-686
【００１２】
【非特許文献４】
M.Koenig etc.:Perfusion CT of the Brain:Radiology
1998;209:85-93
【００１３】
【非特許文献５】
L.Ostergaard etc.:High Resolution Measurement of Cerebral
Blood Flow using Intravascular Tracer Bolus                             Passages:1996;NRM 36:715-725
【００１４】
【発明が解決しようとする課題】
本発明者は、上記従来技術を検討した結果、以下の問題点を見い出した。ガンマ関数解析法は、８〜１０ml／秒程度の高い造影レートが必要であるため侵襲性が高く、定量性も悪いという問題があり、同様に最大勾配法は、８ml／秒程度の高い造影レートが必要であるため侵襲性が高いという問題がある。
【００１５】
一方、deconvolution 法は、３〜５ml／秒程度の低い造影レートで検査可能であるため侵襲性が低く、定量性も優れているが、動脈入力関数と組織残留関数のdeconvolution 演算は非常に演算時間がかかるため、解析結果が出るまでに長時間必要であるという問題がある。特に、急性期脳梗塞のように迅速性を求められる症例においては、血流異常の有無の診断は出来る限り早急に行わなければならないため、演算時間の長さが問題になる。
【００１６】
また、ＣＴ装置を用いて血流動態を検査する場合、パーシャルボリュームアベレージング（以下、ＰＶＡと記述する）効果により動脈を含むボクセルのＣＴ値は本来の値よりも低下する。よって動脈入力関数が本来の値よりも小さくなるため、解析結果の定量性が悪化するという問題がある。
【００１７】
更に、断層像中には、ルームエアーや寝台、骨等の生体機能情報の解析には不要な部分も存在する。断層像中に不要な部分が存在すると、不要領域の信号値の影響で解析精度が悪化するという問題がある。また、不要領域を削除しない場合、不要領域の分だけ無用に演算時間が長くなるという問題もある。更に、不要領域を手動で除去しようとすると、多くの作業時間が消費され、操作者への負担が増大するという問題もある。
【００１８】
本発明の目的は、コンピュータ断層診断装置の提供する断層像から解析対象とする器官の生体機能情報を得る場合において、低い造影レートで定量性のよい生体機能情報を得る事ができ、特に出来る限り短い演算時間で生体機能情報を得る事が可能な血流動態解析装置を提供することにある。
【００１９】
本発明の他の目的は、コンピュータ断層診断装置の提供する断層像から解析対象とする器官の生体機能情報を得る場合において、操作者への負担を増加させることなく解析に不要な領域の信号値の影響による解析精度の悪化を防止する事が可能な血流動態解析装置を提供することにある。
【００２０】
本発明の更に他の目的は、コンピュータ断層診断装置の提供する断層像から解析対象とする器官の生体機能情報を得る場合において、ＰＶＡ効果の影響による定量性の悪化を防止する事が可能な血流動態解析装置を提供することにある。
【００２１】
【課題を解決するための手段】
前記問題を解決するために請求項１に係る血流動態解析装置は、断層像撮影装置により撮影された断層像を入力する入力手段と、この入力手段により入力された断層像の各画素の時間変化の情報を表す時間−濃度曲線を求める第１の演算手段と、この第１の演算手段により求められた断層像の各画素の時間−濃度曲線から流入動脈の時間−濃度曲線を抽出する第２の演算手段と、この第２の演算手段により抽出された流入動脈の時間−濃度曲線から逆フィルタ関数を求める第３の演算手段と、この第３の演算手段により求められた逆フィルタ関数と前記第１の演算手段により求められた断層像の各画素の時間−濃度曲線に基づき前記断層像の各画素の伝達関数を求める第４の演算手段と、この第４の演算手段により求められた前記断層像の各画素の伝達関数を用いて血流動態解析像を求める第５の演算手段と、を備えたことを特徴としている。
【００２２】
即ち、流入動脈の時間−濃度曲線の逆フィルタと各組織における時間−濃度曲線との周波数空間上での掛け算により各組織の伝達関数を求めることができ、従来のdeconvolution 法よりも演算時間の短縮化を図ることができる。
【００２３】
請求項２に係る発明は、前記第３の演算手段は、前記求められた逆フィルタ関数の高周波成分を制御する高周波制御フィルタ手段を備えたことを特徴としている。
また、請求項３に係る発明は、前記高周波制御フィルタ手段は、前記第２の演算手段に抽出された流入動脈の時間―濃度曲線からフィルタ関数を求めるフィルタ関数演算手段を備えたことを特徴としている。
また、請求項４に係る発明は、前記入力手段により入力された断層像のうちの最大連結画素領域を抽出する手段と、この抽出手段により抽出された最大連結画素領域に基づき前記入力手段から入力される断層像のうち血流動態解析に不要な領域を除去する手段と、を備えたことを特徴としている。
また、請求項５に係る発明は、前記不要領域は、ルームエアー、寝台、及び骨を含むことを特徴としている。
また、請求項６に係る発明は、前記入力手段により入力された断層像の各画素の時間−濃度曲線から流出静脈の時間−濃度曲線を求めるとともに、その求められた流入動脈の時間−濃度曲線及び前記流出静脈の時間−濃度曲線のピーク値をそれぞれ求める手段と、前記求められた流入動脈の時間−濃度曲線のピーク値と前記流出静脈の時間―濃度曲線のピーク値とが一致するように流入動脈の時間−濃度曲線におけるパーシャルボリュームアベレージング効果を補正する補正手段と、を備えたことを特徴としている。
【００２４】
また、請求項７に係る血流動態解析方法は、断層像撮影装置により撮影された断層像を入力する入力ステップと、この入力ステップにより入力された断層像の各画素の時間変化の情報を表す時間−濃度曲線を求める第１の演算ステップと、この第１の演算ステップにより求められた断層像の各画素の時間−濃度曲線から流入動脈の時間−濃度曲線を抽出する第２の演算ステップと、この第２の演算ステップにより抽出された流入動脈の時間−濃度曲線から逆フィルタ関数を求める第３の演算ステップと、この第３の演算ステップにより求められた逆フィルタ関数と前記第１の演算ステップにより求められた断層像の各画素の時間−濃度曲線に基づき前記断層像の各画素の伝達関数を求める第４の演算ステップと、この第４の演算ステップにより求められた前記断層像の各画素の伝達関数を用いて血流動態解析像を求める第５の演算ステップと、を含むことを特徴としている。
また、請求項８に係る発明は、前記第３の演算ステップは、前記求められた逆フィルタ関数の高周波成分を制御する高周波制御フィルタステップを含むことを特徴としている。
また、請求項９に係る発明は、前記高周波制御フィルタステップは、前記第２の演算ステップに抽出された流入動脈の時間―濃度曲線からフィルタ関数を求めるフィルタ関数演算ステップを含むことを特徴としている。
また、請求項１０に係る発明は、前記入力ステップにより入力された断層像のうちの最大連結画素領域を抽出するステップと、この抽出ステップにより抽出された最大連結画素領域に基づき前記入力ステップから入力される断層像のうち血流動態解析に不要な領域を除去するステップと、を含むことを特徴としている。
また、請求項１１に係る発明は、前記不要領域は、ルームエアー、寝台、及び骨を含むことを特徴としている。
また、請求項１２に係る発明は、前記入力ステップにより入力された断層像の各画素の時間−濃度曲線から流出静脈の時間−濃度曲線を求めるとともに、その求められた流入動脈の時間−濃度曲線及び前記流出静脈の時間−濃度曲線のピーク値をそれぞれ求めるステップと、前記求められた流入動脈の時間−濃度曲線のピーク値と前記流出静脈の時間―濃度曲線のピーク値とが一致するように流入動脈の時間−濃度曲線におけるパーシャルボリュームアベレージング効果を補正する補正ステップと、を含むことを特徴としている。
また、請求項１３に係る画像診断装置は、生体から断層像を撮影可能な画像診断装置において、前記撮影された断層像を入力する入力手段と、この入力手段により入力された断層像の各画素の時間変化の情報を表す時間−濃度曲線を求める第１の演算手段と、この第１の演算手段により求められた断層像の各画素の時間−濃度曲線から流入動脈の時間−濃度曲線を抽出する第２の演算手段と、この第２の演算手段により抽出された流入動脈の時間−濃度曲線から逆フィルタ関数を求める第３の演算手段と、この第３の演算手段により求められた逆フィルタ関数と前記第１の演算手段により求められた断層像の各画素の時間−濃度曲線に基づき前記断層像の各画素の伝達関数を求める第４の演算手段と、この第４の演算手段により求められた前記断層像の各画素の伝達関数を用いて血流動態解析像を求める第５の演算手段と、を備えたことを特徴としている。
また、請求項１４に係る発明は、前記第３の演算手段は、前記求められた逆フィルタ関数の高周波成分を制御する高周波制御フィルタ手段を備えたことを特徴としている。
また、請求項１５に係る発明は、前記高周波制御フィルタ手段は、前記第２の演算手段に抽出された流入動脈の時間―濃度曲線からフィルタ関数を求めるフィルタ関数演算手段を備えたことを特徴としている。
【００２５】
【発明の実施の形態】
以下、添付図面に従って本発明に係る血流動態解析装置の好ましい実施の形態について詳説する。
【００２６】
図１は本発明に係る血流動態解析装置のハードウェア構成例を示すブロック図である。同図に示すように、この血流動態解析装置は、主としてコンピュータ断層像撮影手段で撮影した断層像を入力する断層像入力装置１と、画像解析などの諸演算を行う演算装置２とから構成されている。
【００２７】
断層像入力装置１と演算装置２は必ずしも別個の装置である必要はなく、一体化していてもよい。また、断層像入力装置１と演算装置２は、必ずしもコンピュータ断層像撮影手段と別個の装置である必要はなく、一体化していてもよい。
【００２８】
演算装置２は、断層像入力装置１から断層像の取り込みを行うインターフェイス（Ｉ／Ｆ）３と、断層像や演算結果を一時的に格納するメモリ４と、諸演算を行う中央処理装置（ＣＰＵ）５と、演算結果や機能画像を記録するハードディスクなどの記録媒体６と、マウスやキーボードなどの外部入力装置７と、処理結果を表示する表示装置８とから構成され、これらは共通バス９によって接続されている。
【００２９】
図２は本発明に係る血流動態解析装置の断層像の入力から機能画像の表示までを示すフローチャートである。
【００３０】
まず、断層像入力装置１から血流動態解析を行う被検体の断層像を選択し、選択された断層像を演算装置２に入力する（ステップ２０１）。入力された断層像はメモリ４に一時格納されるか、または記録媒体６に保存される。断層像入力装置１と演算装置２とが一体化している場合には、外部入力装置７により血流動態解析を行う被検体の断層像を選択し、選択された断層像を記録媒体６から読み出し、メモリ４に一時格納する。
【００３１】
続いて、ＣＰＵ５により、断層像に写り込んでいるルームエアーや寝台、骨等、生体機能情報の解析には不要な領域を取り除く（ステップ２０２）。不要領域除去の方法については後述する。
【００３２】
次に、ＣＰＵ５により、断層像から画素ごとに画素値の時間変化の情報を表す時間−濃度曲線を取り出す（ステップ２０３）。その後、ＣＰＵ５により、時間−濃度曲線から第一循環成分を取り出す（ステップ２０４）。第一循環成分を取り出すアルゴリズムは、ガンマ関数フィッティングや指数関数による外挿等、任意の公知のアルゴリズムでよい。
【００３３】
続いて、ＣＰＵ５により、動脈入力関数に用いる動脈（以下、流入動脈と記述する）とＰＶＡ効果の補正等の定量化に用いる静脈（以下、流出静脈と記述する）を決定する（ステップ２０５）。流入動脈と流出静脈は、操作者が断層像を見ながら手動で指定してもよく、時間−濃度曲線のピーク値及びピーク時間に基づいて解析対象である器官に対する流入動脈と流出静脈を自動的に決定してもよい。
【００３４】
流入動脈と流出静脈を自動的に決定する場合には、まず、流入動脈と流出静脈の中心画素を算出する。流入動脈と流出静脈の中心画素を自動選択するには、各組織における時間−濃度曲線の最大値と最小値の差（以下、ΔTDC と記述する）やピーク時間、ピーク値の特性の違いを利用すればよい。例えば、ΔTDC がある閾値以上の画素の中から、ピーク値が所定の閾値以上で、かつ最もピーク時間の小さい画素を選ぶことにより、流入動脈の中心画素を自動選択することができ、また、ΔTDC がある閾値以上の画素の中からピーク値が所定の閾値以上で、かつ最もピーク時間の大きい画素を選ぶことにより、流出静脈の中心画素を自動選択することができる。次に、流入動脈と流出静脈の中心画素周辺の画素に対して、ピーク値が閾値以上の画素とそうでない画素を分離し、ピーク値が閾値以上の画素のうち、中心画素を含んだ連結画素を抽出することにより、流入動脈領域と流出静脈領域が決定される。尚、中心画素のみを流入動脈と流出静脈としてもよい。
【００３５】
次に、ＣＰＵ５により、ＰＶＡ効果の補正を行う（ステップ２０６）。尚、ＰＶＡ効果の補正方法については後述する。続いて、ＣＰＵ５により、画素毎に逆フィルタ演算を行い、生体機能情報を算出する（ステップ２０７）。尚、逆フィルタ演算については後述する。
【００３６】
画素毎に解析された生体機能情報をマッピングすることで機能画像が得られ（ステップ２０８）、表示装置８により機能画像が表示される（ステップ２０９）。尚、ステップ２０４は、ステップ２０６とステップ２０７の間で実行してもよい。
【００３７】
次に、不要領域を削除する方法について説明する。
【００３８】
図３は不要領域を削除する方法を示すフローチャートである。同図に示すように、まず、ＣＰＵ５により、図２のステップ２０１で入力した断層像を２値化する（ステップ３０１）。２値化処理は、例えば、断層像において、閾値以上の画素値を有する画素を１に、閾値未満の画素値を有する画素を０にそれぞれ置き換えることで達成される。尚、閾値は、解析対象となる生体組織とルームエアーが分離できるような任意の値でよい。例えば、入力した断層像がＣＴ像であるならば、閾値を−２００程度に設定すればよい。
【００３９】
続いて、ＣＰＵ５により、ステップ３０１で２値化された画像をラベリング処理する（ステップ３０２）。ラベリング処理は、つながっている全ての画素（連結成分）に同じ番号（ラベル）を付けて、異なる連結成分には別の番号を付けることで達成される。尚、隣接する画素が互いに１の場合には、これらの画素はつながっていることになる。ラベルとして付与する値は、例えば５０からの続き番号とする。
【００４０】
次に、ＣＰＵ５により、最大連結成分を探索する（ステップ３０３）。最大連結成分の探索は、ラベリングした画像全体を走査して、ラベル値ごとに画素数を数える処理を行い、最も画素数が多いラベル値を有する連結成分を選び出すことで達成される。その後、ＣＰＵ５により、最大連結成分を抽出する（ステップ３０４）。最大連結成分の抽出は、ステップ３０３で選び出された最も画素数が多いラベル値を有する画素を残し、それ以外の画素を０に置き換えることで達成される。これにより２値化処理の段階では画素値が１に置き換えられる可能性のある寝台や点滴チューブ等の不要領域の画素値が０に置き換えられる。
【００４１】
次に、ＣＰＵ５により、最大連結成分の輪郭を追跡し、最大連結成分の最外周の輪郭線を抽出する（ステップ３０５）。輪郭線の追跡は、ステップ３０４で最大連結成分のみが残された画像の左上隅画素から横に走査し、最初に出会った非ゼロのラベル値を有する画素を開始点として反時計回り方向に輪郭を追跡し、開始点に戻ったら追跡を終了することで達成される。上記輪郭線上の画素にはラベル値とは異なる値、例えば１を代入しておく。
【００４２】
続いて、ＣＰＵ５により、ステップ３０５で求めた輪郭線の内部を１で塗りつぶし、マスク画像を作成する（ステップ３０６）。輪郭線の内部を塗りつぶすには、従来から知られている、閉領域内部の１点を開始点として閉領域内部を塗りつぶす処理であるシードフィルアルゴリズムを適用すればよい。本実施の形態では、閉領域はステップ３０５で求めた輪郭線である。また、閉領域内部の点は、ステップ３０４で求めた最大連結成分のラベル値の付いた画素である。このラベル値の付いた画素を１つ検出し、その点を開始点としてシードフィル処理を行えば、輪郭線内部を塗りつぶすことができる。
【００４３】
次に、ＣＰＵ５により、入力画像（断層像）とステップ３０６で作成したマスク画像とを掛け合わせ、マスク領域外の画素を削除する（ステップ３０７）。具体的には、マスク画像を１画素ずつ調べ、ゼロの画素を検出したら、その画素の座標に対応する原画像の画素を最低輝度に置き換える処理をする。これを全画素について行えば、マスク領域外の画素が入力画像から削除される。
【００４４】
次に、ＰＶＡ効果の補正方法について説明する。
【００４５】
図４は本実施の形態におけるＰＶＡ効果の補正方法を示す図である。例えば、脳の生体機能情報を解析する場合、一般に流入動脈には前大脳動脈や中大脳動脈を、流出静脈には上矢状静脈洞を用いる。Lapin らの研究(J Comput Assisit Tomogr 1993;108-114)によれば、ＰＶＡ効果により動脈を含むボクセルの信号値が低下するため、流入動脈である前大脳動脈や中大脳動脈の時間−濃度曲線は本来の値よりも小さくなる。またLee らの研究(AJNR 1999;20:63-73)によれば、直径1.73mm以上の血管はＰＶＡ効果の影響を受けないことから、流出静脈である上矢状静脈洞の時間−濃度曲線はＰＶＡ効果の影響を受けていない。
【００４６】
よって流入動脈と流出静脈の時間濃度曲線は、一般に図４（ａ）のようになり、流入動脈の時間濃度曲線のピーク値Ｐａは、流出静脈のピーク値ＰｖよりもＰＶＡ効果による信号値の低下の分だけ小さな値になる。理想的には血管内の造影剤濃度は血管径の大小に依らず一定であるので、ＰＶＡ効果の影響が無ければ信号値のピーク値は血管径の大小によらず一定である。
【００４７】
従って、流入動脈のピーク値をＰａ、流出静脈のピーク値をＰｖとすると、図４（ｂ）に示すように流入動脈の時間−濃度曲線をＰｖ／Ｐａ倍すれば、流入動脈におけるＰＶＡ効果を補正できる。ＰＶＡ効果を補正する前の時間ｔにおける流入動脈の時間−濃度曲線をAIFpva(t) 、ＰＶＡ効果を補正した彼の時間ｔにおける流入動脈の時間−濃度曲線をAIF(t)とすると、両者の間に次式が成り立つように補正すればよい。
【００４８】
【数１】
AIF(t)＝AIFpva(t) ×Pv／Pa      …（１）
図５は本実施の形態におけるＰＶＡ効果の他の補正方法を示す図である。
【００４９】
過去の臨床画像の中に、これから実施する検査と全く同一の撮影条件、造影条件で撮影した例があるならば、この過去の臨床画像を用いてＰＶＡ効果を補正することも可能である。
【００５０】
例えば、図５（ａ）に示す臨床画像において、断層像内に含まれる複数の血管１〜５に対し、血管径Ｒ１〜Ｒ５とその血管の時間−濃度曲線のピーク値を算出する。これを図表に表すと、図５（ｂ）となり、プロットしてグラフに表すと、図４（ｃ）のようなカーブが得られる。尚、血管径とピーク値の算出方法については後述する。血管径が細いほどＰＶＡ効果の影響を強く受け、血管径がある値以上であればＰＶＡ効果の影響を受けないため、図４（ｃ）のカーブはある径Ｒthまでは徐々に上昇し、Ｒth以上のＲに対してはピーク値が平衡になる。尚、この平衡状態におけるピーク値を便宜上Ｐthとし、流入動脈における血管径をＲａ、流入動脈における時間−濃度曲線のピーク値をＰａとする。
【００５１】
ＲａがＲthより小さい場合には、流入動脈の時間−濃度曲線をＰth／Ｐａ倍すれば、流入動脈におけるＰＶＡ効果を補正できる。ＰＶＡ効果を補正する前の時間ｔにおける流入動脈の時間−濃度曲線をAIFpva(t) 、ＰＶＡ効果を補正した後の時間ｔにおける流入動脈の時間−濃度曲線をAIF(t)とすると、両者の間に次式が成り立つように補正すればよい。
【００５２】
【数２】
AIF(t)＝AIFpva(t) ×Pth／Pa      …（２）
次に、血管径とピーク値の算出方法について説明する。
【００５３】
図６は本実施の形態における血管径とピーク値の算出手順を示すフローチャートである。まず、ＣＰＵ５により、各時間−濃度曲線のピーク時間とピーク値を算出する（ステップ６０１）。次に、ＣＰＵ５により、血管の中心画素を算出する（ステップ６０２）。中心画素を算出するには、ピーク値の特性を利用すればよく、ピーク値が最大の画素を選ぶことにより、血管の中心画素を自動選択することができる。続いて、ＣＰＵ５により、血管の中心画素周辺の画素に対して、ピーク値が閾値以上の画素とそうでない画素を分離する（ステップ６０３）。その後、ステップ６０３で求めたピーク値が閾値以上の画素のうち、中心画素を含んだ連結画素を抽出することにより、血管領域が決定される（ステップ６０４）。
【００５４】
次に、ＣＰＵ５により、血管領域内の画素数を算出し（ステップ６０５）、その算出した画素数から血管領域の面積を算出する（ステップ６０６）。一画素の面積は装置毎に一定値であるので、一画素あたりの面積にステップ６０５で求められた画素数を掛けることで血管領域の面積Ｓが得られる。
【００５５】
続いて、ＣＰＵ５により、血管径Ｒを算出する（ステップ６０７）。断層像上での血管の形状は円形で近似して差し支えないので、血管径Ｒと面積Ｓには、次式、
【００５６】
【数３】
Ｒ＝√（Ｓ／π）      …（３）
が成り立ち、（３）式に従って血管径を算出すればよい。
【００５７】
次に、ＣＰＵ５により血管の時間−濃度曲線を算出する（ステップ６０８）。ステップ２０３において既に各画素における時間−濃度曲線が算出されているので、血管領域の全ての画素における時間−濃度曲線の平均値を算出することで、血管の時間−濃度曲線が得られる。また、平均値をとるかわりに、ステップ６０２で算出されている血管中心画素からの距離に応じて、任意に重み付けをしてもよい。
【００５８】
次に、ＣＰＵ５により、血管の時間−濃度曲線のピーク値を算出する（ステップ６０９）。尚、ステップ６０５からステップ６０７までの処理と、ステップ６０８からステップ６０９までの処理は、処理の順序が逆になってもよい。本実施例では、ピーク値の特性を利用して血管の抽出を行ったが、血管の抽出法はこれに限定されるものではなく、血管が正確に抽出できる方法であるならば任意の方法でよい。また、断層像を見ながら手動で血管領域を指定することにより、血管を抽出してもよい。
【００５９】
次に、逆フィルタ法の概念について説明する。
【００６０】
図７は本実施の形態における逆フィルタ法を説明するための概念図である。伝達関数ｆの対象組織において、入力関数をCin(t)、出力関数をCout(t) とおくと、Cout(t) はCin(t)とｆを畳み込み積分したものになり、次式、
【００６１】
【数４】
0004013128.tif 000002
のように表される。
【００６２】
造影剤を動脈から体内に注入する場合、Cin(t)は流入動脈における時間−濃度曲線（動脈入力関数）AIF(t)に、Cout(t) は組織における時間−濃度曲線Cr(t) に置き換えられ、次式、
【００６３】
【数５】
0004013128.tif 000003
のように表される。
【００６４】
（５）式を周波数空間に変換すると、次式、
【００６５】
【数６】
Ｆ｛Cr(t)｝＝Ｆ｛AIF(t)｝・Ｆ｛f(t)｝…（６）
のようになる。
【００６６】
ここで、Ｆ｛AIF(t)｝はAIF(t)のフーリエ変換を、Ｆ｛Cr(t) ｝はCr(t) のフーリエ変換を、Ｆ｛f(t)｝はf(t)のフーリエ変換をそれぞれ表す。
【００６７】
（５）式と（６）式より、流入動脈における時間−濃度曲線と組織における時間−濃度曲線を逆畳み込み演算すれば伝達関数ｆを求めることができるが、逆畳み込み演算を行う代わりに、流入動脈における時間−濃度曲線の逆フィルタを組織の時間−濃度曲線に掛けても等価な演算結果が得られる。
【００６８】
AIF(t)の逆フィルタをAIF(t)-1とおくと、（６）式より、f(t)は次式、
【００６９】
【数７】
f(t)＝Ｆ-1［Ｆ｛AIF-1(t) ｝・Ｆ｛Cr(t) ｝］        …（７）
のように表される。
【００７０】
ここでＦ-1は、逆フーリエ変換を表す。
【００７１】
図８は本実施の形態における逆フィルタ法を説明するための概念図である。
【００７２】
血流動態解析を行う場合、被検者への侵襲性を抑えるため、たいていの場合、造影剤は静脈から体内へ注入する。この場合、造影剤は体内を循環した後、解析対象組織に流入する。入力関数をCin(t)、循環経路における伝達関数をｆ1 、対象組織における伝達関数をｆ2 、対象組織からの出力関数をCout(t) とおくと、Cout(t) は、Cin(t)とｆ1 を畳み込み積分したものを、更にｆ2 と畳み込み積分したものになるので、次式、
【００７３】
【数８】
0004013128.tif 000004
のようになる。
【００７４】
Cin(t)とｆ1 のたたみこみ積分は流入動脈における時間−濃度曲線（動脈入力関数）AIF(t)に、Cout(t) は組織における時間−濃度曲線Cr(t) に置き換えられ、次式、
【００７５】
【数９】
0004013128.tif 000005
のように表される。
【００７６】
（９）式は、前述した（５）式と等価である。よって伝達関数ｆ2 は、次式、
【００７７】
【数１０】
f2(t) ＝Ｆ-1［Ｆ｛AIF-1(t) ｝・Ｆ｛Cr(t) ｝］    …（１０）
のように表される。
【００７８】
次に、逆フィルタ演算について説明する。
【００７９】
図９は本実施の形態における逆フィルタ演算の処理手順を示すフローチャートである。同図に示すように、ＣＰＵ５により、流入動脈の時間−濃度曲線をフーリエ変換する（ステップ９０１）。次に、ＣＰＵ５により、流入動脈の時間−濃度曲線のフーリエ変換から逆フィルタを算出する（ステップ９０２）。流入動脈の時間−濃度曲線のフーリエ変換の実部をＦ｛AIFR｝、虚部をＦ｛AIFI｝とすると、逆フィルタは、次式、
【００８０】
【数１１】
0004013128.tif 000006
で表される。
【００８１】
次に、逆フィルタの高周波数成分を制御する（ステップ９０３）。ウィナーフィルタやバターワースフィルタなどの高周波数領域制御フィルタによって逆フィルタをフィルタリングすることにより、逆フィルタの高周波数成分が制御される。
【００８２】
続いて、ＣＰＵ５により、各組織の時間−濃度曲線をフーリエ変換する（ステップ９０４）。このフーリエ変換に対して、ステップ９０２で求めた逆フィルタでフィルタリングする（ステップ９０５）。このフィルタリングされた組織の時間−濃度曲線のフーリエ変換を逆フーリエ変換する伝達関数が求められる（ステップ９０６）。
【００８３】
次に、ＣＰＵ５により、ステップ９０６で求められた伝達関数から、生体機能情報を算出する（ステップ９０７）。生体機能情報の算出方法については後述する。尚、ステップ９０４からステップ９０７までは、解析する組織の断層像上での画素数分だけ繰り返す。また、ステップ９０１からステップ９０３までと、ステップ９０４は、処理の順序は逆になってもよい。
【００８４】
次に、伝達関数から生体機能情報を算出する方法について説明する。
【００８５】
図１０は生体機能情報の算出方法の説明図である。伝達関数をf(t)とすると、f(t)は、例えば図１０のような形状になり、血流量（bloodflow:以下、BFと記述する）は、次式、
【００８６】
【数１２】
BF＝fmax         …（１２）
のようにf(t)の最大値から求められ、血液量（bloodvolume:以下、BVと記述する）は、次式、
【００８７】
【数１３】
BV＝∫f(t)dt      …（１３）
のようにf(t)の曲線下面積から求められ、平均通過時間（meantransit time: 以下、MTT と記述する）は、次式、
【００８８】
【数１４】
MTT ＝∫f(t)dt／fmax     …（１４）
のようにf(t)の幅から求められる。
【００８９】
ここで、f max は、f(t)の最大値を表す。
【００９０】
【発明の効果】
以上説明したように本発明によれば、コンピュータ断層診断装置の提供するダイナミック撮影した断層像から生体機能情報解析を行う場合、ガンマ関数解析法や最大勾配法よりも低い造影レートで、deconvolution 法よりも短い演算時間で、生体機能情報解析を行えるという効果がある。また、ＰＶＡ効果による定量性の悪化を防止できるという効果もある。更に、操作者への負担を増加させることなく不要領域の影響による無用な演算時間の増加や解析精度の悪化を防止する事ができるという効果もある。
【図面の簡単な説明】
【図１】  本発明に係る血流動態解析装置のハードウェア構成例を示すブロック図
【図２】  本発明に係る血流動態解析装置の断層像の入力から機能画像の表示までを示すフローチャート
【図３】  不要領域を削除する方法を示すフローチャート
【図４】  本実施の形態におけるＰＶＡ効果の補正方法を示す図
【図５】  本実施の形態におけるＰＶＡ効果の他の補正方法を示す図
【図６】  本実施の形態における血管径とピーク値の算出手順を示すフローチャート
【図７】  本実施の形態における逆フィルタ法を説明するための概念図
【図８】  本実施の形態における逆フィルタ法を説明するための概念図
【図９】  本実施の形態における逆フィルタ演算の処理手順を示すフローチャート
【図１０】  生体機能情報の算出方法の説明図
【符号の説明】
１…断層像入力装置、２…演算装置、３…インターフェイス（Ｉ／Ｆ）、４…モリ、５…中央処理装置（ＣＰＵ）、６…記録媒体、７…外部入力装置、８…表示装置、９…共通バス
【図１】
0004013128.tif 000007
【図２】
0004013128.tif 000008
【図３】
0004013128.tif 000009
【図４】
0004013128.tif 000010
【図５】
0004013128.tif 000011
【図６】
0004013128.tif 000012
【図７】
0004013128.tif 000013
【図８】
0004013128.tif 000014
【図９】
0004013128.tif 000015
【図１０】
0004013128.tif 000016