JP2004097665A - 血流動態解析装置 - Google Patents
血流動態解析装置 Download PDFInfo
- Publication number
- JP2004097665A JP2004097665A JP2002266475A JP2002266475A JP2004097665A JP 2004097665 A JP2004097665 A JP 2004097665A JP 2002266475 A JP2002266475 A JP 2002266475A JP 2002266475 A JP2002266475 A JP 2002266475A JP 2004097665 A JP2004097665 A JP 2004097665A
- Authority
- JP
- Japan
- Prior art keywords
- time
- curve
- blood flow
- tomographic image
- concentration curve
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000004458 analytical method Methods 0.000 title claims abstract description 33
- 230000017531 blood circulation Effects 0.000 title claims abstract description 29
- 210000001367 artery Anatomy 0.000 claims abstract description 48
- 230000006870 function Effects 0.000 claims abstract description 30
- 230000008827 biological function Effects 0.000 claims abstract description 29
- 238000012546 transfer Methods 0.000 claims abstract description 17
- 230000000694 effects Effects 0.000 claims description 32
- 210000001519 tissue Anatomy 0.000 claims description 28
- 210000003462 vein Anatomy 0.000 claims description 23
- 238000012937 correction Methods 0.000 claims description 10
- 210000000988 bone and bone Anatomy 0.000 claims description 4
- 238000003384 imaging method Methods 0.000 abstract description 11
- 210000000056 organ Anatomy 0.000 abstract description 8
- 238000003745 diagnosis Methods 0.000 abstract description 3
- 230000009466 transformation Effects 0.000 abstract 3
- 238000000034 method Methods 0.000 description 48
- 210000004204 blood vessel Anatomy 0.000 description 35
- 238000004364 calculation method Methods 0.000 description 13
- 238000010586 diagram Methods 0.000 description 12
- 238000002591 computed tomography Methods 0.000 description 10
- 238000012905 input function Methods 0.000 description 9
- 239000002872 contrast media Substances 0.000 description 8
- 238000012545 processing Methods 0.000 description 6
- 239000000284 extract Substances 0.000 description 5
- 230000006866 deterioration Effects 0.000 description 4
- 210000004556 brain Anatomy 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 238000012935 Averaging Methods 0.000 description 2
- 230000001154 acute effect Effects 0.000 description 2
- 210000002551 anterior cerebral artery Anatomy 0.000 description 2
- 239000008280 blood Substances 0.000 description 2
- 210000004369 blood Anatomy 0.000 description 2
- 230000003727 cerebral blood flow Effects 0.000 description 2
- 206010008118 cerebral infarction Diseases 0.000 description 2
- 208000026106 cerebrovascular disease Diseases 0.000 description 2
- 230000004087 circulation Effects 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 230000000004 hemodynamic effect Effects 0.000 description 2
- 238000010191 image analysis Methods 0.000 description 2
- 210000003657 middle cerebral artery Anatomy 0.000 description 2
- 238000009206 nuclear medicine Methods 0.000 description 2
- 230000010412 perfusion Effects 0.000 description 2
- 238000011002 quantification Methods 0.000 description 2
- 210000000798 superior sagittal sinus Anatomy 0.000 description 2
- ZCYVEMRRCGMTRW-UHFFFAOYSA-N 7553-56-2 Chemical compound [I] ZCYVEMRRCGMTRW-UHFFFAOYSA-N 0.000 description 1
- 229910052688 Gadolinium Inorganic materials 0.000 description 1
- 208000006011 Stroke Diseases 0.000 description 1
- 230000002159 abnormal effect Effects 0.000 description 1
- 235000013339 cereals Nutrition 0.000 description 1
- 238000002405 diagnostic procedure Methods 0.000 description 1
- 238000002597 diffusion-weighted imaging Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 238000013213 extrapolation Methods 0.000 description 1
- UIWYJDYFSGRHKR-UHFFFAOYSA-N gadolinium atom Chemical compound [Gd] UIWYJDYFSGRHKR-UHFFFAOYSA-N 0.000 description 1
- 229910052740 iodine Inorganic materials 0.000 description 1
- 239000011630 iodine Substances 0.000 description 1
- 238000002372 labelling Methods 0.000 description 1
- 230000003902 lesion Effects 0.000 description 1
- 230000005291 magnetic effect Effects 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 239000003973 paint Substances 0.000 description 1
- 238000002600 positron emission tomography Methods 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 239000000700 radioactive tracer Substances 0.000 description 1
- 238000002603 single-photon emission computed tomography Methods 0.000 description 1
- 238000003325 tomography Methods 0.000 description 1
Images
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/02—Detecting, measuring or recording pulse, heart rate, blood pressure or blood flow; Combined pulse/heart-rate/blood pressure determination; Evaluating a cardiovascular condition not otherwise provided for, e.g. using combinations of techniques provided for in this group with electrocardiography or electroauscultation; Heart catheters for measuring blood pressure
- A61B5/026—Measuring blood flow
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/48—Diagnostic techniques
- A61B6/481—Diagnostic techniques involving the use of contrast agents
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/50—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications
- A61B6/504—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications for diagnosis of blood vessels, e.g. by angiography
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/50—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications
- A61B6/507—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications for determination of haemodynamic parameters, e.g. perfusion CT
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
- G06T7/0012—Biomedical image inspection
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7235—Details of waveform analysis
- A61B5/7253—Details of waveform analysis characterised by using transforms
- A61B5/7257—Details of waveform analysis characterised by using transforms using Fourier transforms
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30101—Blood vessel; Artery; Vein; Vascular
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Medical Informatics (AREA)
- General Health & Medical Sciences (AREA)
- Physics & Mathematics (AREA)
- Public Health (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Biophysics (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Radiology & Medical Imaging (AREA)
- Molecular Biology (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- Pathology (AREA)
- Veterinary Medicine (AREA)
- Optics & Photonics (AREA)
- High Energy & Nuclear Physics (AREA)
- Dentistry (AREA)
- Oral & Maxillofacial Surgery (AREA)
- Cardiology (AREA)
- Physiology (AREA)
- Vascular Medicine (AREA)
- Hematology (AREA)
- Quality & Reliability (AREA)
- Computer Vision & Pattern Recognition (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Apparatus For Radiation Diagnosis (AREA)
- Image Processing (AREA)
- Magnetic Resonance Imaging Apparatus (AREA)
Abstract
【課題】コンピュータ断層診断装置の提供する断層像から解析対象とする器官の生体機能情報を得る場合において、低い造影レートで定量性のよい生体機能情報を得る事ができ、特に出来る限り短い演算時間で生体機能情報を得る事が可能な血流動態解析装置を提供する。
【解決手段】ダイナミック撮影した断層像に基づいて流入動脈の時間−濃度曲線及び各組織における時間−濃度曲線をフーリエ変換する（ステップ９０１、９０４）。流入動脈の時間−濃度曲線のフーリエ変換から逆フィルタを算出する（ステップ９０２）。この逆フィルタを各組織における時間−濃度曲線のフーリエ変換に掛け合わせることで各組織の伝達関数を求める（ステップ９０５、９０６）。このようにして求めた各組織の伝達関数を使用して生体機能情報を算出する（ステップ９０７）。
【選択図】 図９
【解決手段】ダイナミック撮影した断層像に基づいて流入動脈の時間−濃度曲線及び各組織における時間−濃度曲線をフーリエ変換する（ステップ９０１、９０４）。流入動脈の時間−濃度曲線のフーリエ変換から逆フィルタを算出する（ステップ９０２）。この逆フィルタを各組織における時間−濃度曲線のフーリエ変換に掛け合わせることで各組織の伝達関数を求める（ステップ９０５、９０６）。このようにして求めた各組織の伝達関数を使用して生体機能情報を算出する（ステップ９０７）。
【選択図】 図９
Description
【０００１】
【発明の属する技術分野】
本発明は血流動態解析装置に係り、特にＸ線や電磁波を用いて被検体の時間変化に依存する生体機能情報を得ることが可能なコンピュータ断層診断装置の提供する断層像から生体機能情報解析等の画像解析を行うことに対して有用な技術に関する。
【０００２】
【従来の技術】
従来、ダイナミック撮影を行うコンピュータ断層診断装置として、陽電子放射断層撮影装置（Ｐｏｓｉｔｒｏｎ Ｅｍｉｓｓｉｏｎ Ｔｏｍｏｇｒａｐｈｙ：ＰＥＴ）や単光子放射ＣＴ（ＳｉｎｇｌｅＰｈｏｔｏｎ Ｅｍｉｓｓｉｏｎ ＣＴ：ＳＰＥＣＴ）などの核医学装置がある。核医学装置による血流動態計測は、放射性核種を被検体内に注入し、核種からの放射線をシンチレーションカメラで計測すると、横断面の放射性核種分布像を得られ、この放射性核種分布像を解析することで、器官の生体機能情報を解析していた。
【０００３】
他のダイナミック撮影を行うコンピュータ断層診断装置として、Ｘ線ＣＴ装置がある。単純Ｘ線ＣＴ像では病変の超早期相の診断が困難である場合には、Ｘ線ＣＴ装置によるダイナミック撮影を行い、この撮影時にコントラストを付けるために主にヨード系の造影剤を用いていた。
【０００４】
造影剤を注入してダイナミックスキャンを行うと、造影剤濃度とＣＴ値は比例関係にあることから、撮影断面の時間変化の情報が得られる。撮影断面の時間変化から各組織のＣＴ値の時間変化である時間−濃度曲線が得られ、この時間−濃度曲線を解析することで器官の生体機能情報を解析していた（例えば非特許文献１参照。）。
【０００５】
更に他のダイナミック撮影を行うコンピュータ断層診断装置として、ＭＲ装置がある。高磁場ＭＲ装置では非造影の診断方法もあるが、高磁場ＭＲ装置でコントラストの良い画像を得たい場合や低磁場ＭＲ装置においては、主にガドリニウム系の造影剤を用いていた。ＭＲ装置においても撮影断面の時間変化から各組織のＭＲ値の相対的信号強度変化率である時間−濃度曲線が得られる。この時間−濃度曲線を解析することで器官の時間変化に依存する生体機能情報を解析していた（例えば非特許文献２参照。）。
【０００６】
一方、血流動態を解析するアルゴリズムとして、ガンマ関数解析法がある。ガンマ関数解析法では、時間−濃度曲線をガンマ関数でフィッティングし、フィッティングカーブのピーク値や曲線下面積から血流情報を算出していた（例えば非特許文献３参照。）。
【０００７】
他の血流動態を解析するアルゴリズムとして、最大勾配法（ｍａｘｉｍｕｍ ｓｌｏｐｅ 法）がある。最大勾配法では、各組織における時間−濃度曲線の傾きの最大値を動脈入力関数におけるＣＴ値上昇の最大値で割ることで血流量を算出していた（例えば非特許文献４参照。）。
【０００８】
更に他の血流動態を解析するアルゴリズムとして、ｄｅｃｏｎｖｏｌｕｔｉｏｎ（デコンボリューション）法がある。ｄｅｃｏｎｖｏｌｕｔｉｏｎ 法では、動脈入力関数と組織残留関数をｄｅｃｏｎｖｏｌｕｔｉｏｎ 演算することでインパルス残留関数を求め、インパルス残留関数のピーク値や曲線下面積から血流情報を算出していた（例えば非特許文献５参照。）。
【０００９】
【非特許文献１】
Ｊ．Ｐ．Ｍｕｉｚｅｌａａｒ ｅｔｃ．：Ａ Ｎｅｗ Ｍｅｔｈｏｄ ｆｏｒ Ｑｕａｎｔｉｔａｔｉｖｅ Ｒｅｇｉｏｎａｌ Ｃｅｒｅｂｒａｌ ＢｌｏｏｄＶｏｌｕｍｅ Ｍｅａｓｕｒｅｍｅｎｔｓ Ｕｓｉｎｇ Ｃｏｍｐｕｔｅｄ Ｔｏｍｏｇｒａｐｈｙ．：Ｓｔｒｏｋｅ１９９７；２８：１９９８−２００５
【００１０】
【非特許文献２】
井川房雄 他：拡散強調画像と潅流画像による急性期脳梗塞の診断の意義：日磁医誌 第２０巻５号２０３−２０９頁（２０００年）
【００１１】
【非特許文献３】
Ｌ．Ａｘｅｌ：Ｃｅｒｅｂｒａｌ Ｂｌｏｏｄ Ｆｌｏｗ Ｄｅｔｅｒｍｉｎａｔｉｏｎ ｂｙ Ｒａｐｉｄ−Ｓｅｑｕｅｎｃｅ Ｃｏｍｐｕｔｅｄ Ｔｏｍｏｇｒａｐｈｙ：Ｒａｄｉｏｌｏｇｙ １９８０；１３７：６７９−６８６
【００１２】
【非特許文献４】
Ｍ．Ｋｏｅｎｉｇ ｅｔｃ．：Ｐｅｒｆｕｓｉｏｎ ＣＴ ｏｆ ｔｈｅ Ｂｒａｉｎ：Ｒａｄｉｏｌｏｇｙ
１９９８；２０９：８５−９３
【００１３】
【非特許文献５】
Ｌ．Ｏｓｔｅｒｇａａｒｄ ｅｔｃ．：Ｈｉｇｈ Ｒｅｓｏｌｕｔｉｏｎ Ｍｅａｓｕｒｅｍｅｎｔ ｏｆ Ｃｅｒｅｂｒａｌ Ｂｌｏｏｄ Ｆｌｏｗ ｕｓｉｎｇ Ｉｎｔｒａｖａｓｃｕｌａｒ Ｔｒａｃｅｒ Ｂｏｌｕｓ Ｐａｓｓａｇｅｓ：１９９６；ＮＲＭ ３６：７１５−７２５
【００１４】
【発明が解決しようとする課題】
本発明者は、上記従来技術を検討した結果、以下の問題点を見い出した。ガンマ関数解析法は、８〜１０ｍｌ／秒程度の高い造影レートが必要であるため侵襲性が高く、定量性も悪いという問題があり、同様に最大勾配法は、８ｍｌ／秒程度の高い造影レートが必要であるため侵襲性が高いという問題がある。
【００１５】
一方、ｄｅｃｏｎｖｏｌｕｔｉｏｎ 法は、３〜５ｍｌ／秒程度の低い造影レートで検査可能であるため侵襲性が低く、定量性も優れているが、動脈入力関数と組織残留関数のｄｅｃｏｎｖｏｌｕｔｉｏｎ 演算は非常に演算時間がかかるため、解析結果が出るまでに長時間必要であるという問題がある。特に、急性期脳梗塞のように迅速性を求められる症例においては、血流異常の有無の診断は出来る限り早急に行わなければならないため、演算時間の長さが問題になる。
【００１６】
また、ＣＴ装置を用いて血流動態を検査する場合、パーシャルボリュームアベレージング（以下、ＰＶＡと記述する）効果により動脈を含むボクセルのＣＴ値は本来の値よりも低下する。よって動脈入力関数が本来の値よりも小さくなるため、解析結果の定量性が悪化するという問題がある。
【００１７】
更に、断層像中には、ルームエアーや寝台、骨等の生体機能情報の解析には不要な部分も存在する。断層像中に不要な部分が存在すると、不要領域の信号値の影響で解析精度が悪化するという問題がある。また、不要領域を削除しない場合、不要領域の分だけ無用に演算時間が長くなるという問題もある。更に、不要領域を手動で除去しようとすると、多くの作業時間が消費され、操作者への負担が増大するという問題もある。
【００１８】
本発明の目的は、コンピュータ断層診断装置の提供する断層像から解析対象とする器官の生体機能情報を得る場合において、低い造影レートで定量性のよい生体機能情報を得る事ができ、特に出来る限り短い演算時間で生体機能情報を得る事が可能な血流動態解析装置を提供することにある。
【００１９】
本発明の他の目的は、コンピュータ断層診断装置の提供する断層像から解析対象とする器官の生体機能情報を得る場合において、操作者への負担を増加させることなく解析に不要な領域の信号値の影響による解析精度の悪化を防止する事が可能な血流動態解析装置を提供することにある。
【００２０】
本発明の更に他の目的は、コンピュータ断層診断装置の提供する断層像から解析対象とする器官の生体機能情報を得る場合において、ＰＶＡ効果の影響による定量性の悪化を防止する事が可能な血流動態解析装置を提供することにある。
【００２１】
【課題を解決するための手段】
上記問題を解決するために請求項１に係る血流動態解析装置は、コンピュータ断層像撮影手段でダイナミック撮影した断層像の入力手段と、前記入力手段から入力した断層像に基づいて各画素ごとに画素値の時間変化の情報を表す時間−濃度曲線を求める演算手段であって、流入動脈の時間−濃度曲線及び各組織における時間−濃度曲線を求める演算手段と、前記求めた流入動脈の時間−濃度曲線の逆フィルタを算出する手段と、前記算出した逆フィルタを前記求めた各組織における時間−濃度曲線に掛け合わせることで各組織の伝達関数を求める手段と、各組織の伝達関数を使用して生体機能情報を算出する手段と、を備えたことを特徴としている。
【００２２】
即ち、流入動脈の時間−濃度曲線の逆フィルタと各組織における時間−濃度曲線との周波数空間上での掛け算により各組織の伝達関数を求めることができ、従来のｄｅｃｏｎｖｏｌｕｔｉｏｎ 法よりも演算時間の短縮化を図ることができる。
【００２３】
請求項２に示すように前記演算手段は、断層像内の最大連結画素領域を抽出する手段と、前記抽出した最大連結画素領域を示すマスク画像を作成する手段と、前記マスク画像を使用して前記入力手段から入力する断層像のうちの当該マスク画像の領域以外の画素を削除する手段とを有し、ルームエアーや寝台、骨等の生体機能情報解析には不要な領域を自動的に除去することを特徴としている。
【００２４】
請求項３に示すように前記演算手段は、流出静脈の時間−濃度曲線を求めるとともに、前記求めた流入動脈の時間−濃度曲線及び流出静脈の時間−濃度曲線のピーク値をそれぞれ求める手段と、前記流入動脈の時間−濃度曲線のピーク値が前記流出静脈の時間−濃度曲線のピーク値が一致するように前記流入動脈の時間−濃度曲線を補正する補正手段とを有し、流入動脈の時間−濃度曲線におけるパーシャルボリュームアベレージング効果を補正することを特徴としている。
【００２５】
【発明の実施の形態】
以下、添付図面に従って本発明に係る血流動態解析装置の好ましい実施の形態について詳説する。
【００２６】
図１は本発明に係る血流動態解析装置のハードウェア構成例を示すブロック図である。同図に示すように、この血流動態解析装置は、主としてコンピュータ断層像撮影手段で撮影した断層像を入力する断層像入力装置１と、画像解析などの諸演算を行う演算装置２とから構成されている。
【００２７】
断層像入力装置１と演算装置２は必ずしも別個の装置である必要はなく、一体化していてもよい。また、断層像入力装置１と演算装置２は、必ずしもコンピュータ断層像撮影手段と別個の装置である必要はなく、一体化していてもよい。
【００２８】
演算装置２は、断層像入力装置１から断層像の取り込みを行うインターフェイス（Ｉ／Ｆ）３と、断層像や演算結果を一時的に格納するメモリ４と、諸演算を行う中央処理装置（ＣＰＵ）５と、演算結果や機能画像を記録するハードディスクなどの記録媒体６と、マウスやキーボードなどの外部入力装置７と、処理結果を表示する表示装置８とから構成され、これらは共通バス９によって接続されている。
【００２９】
図２は本発明に係る血流動態解析装置の断層像の入力から機能画像の表示までを示すフローチャートである。
【００３０】
まず、断層像入力装置１から血流動態解析を行う被検体の断層像を選択し、選択された断層像を演算装置２に入力する（ステップ２０１）。入力された断層像はメモリ４に一時格納されるか、または記録媒体６に保存される。断層像入力装置１と演算装置２とが一体化している場合には、外部入力装置７により血流動態解析を行う被検体の断層像を選択し、選択された断層像を記録媒体６から読み出し、メモリ４に一時格納する。
【００３１】
続いて、ＣＰＵ５により、断層像に写り込んでいるルームエアーや寝台、骨等、生体機能情報の解析には不要な領域を取り除く（ステップ２０２）。不要領域除去の方法については後述する。
【００３２】
次に、ＣＰＵ５により、断層像から画素ごとに画素値の時間変化の情報を表す時間−濃度曲線を取り出す（ステップ２０３）。その後、ＣＰＵ５により、時間−濃度曲線から第一循環成分を取り出す（ステップ２０４）。第一循環成分を取り出すアルゴリズムは、ガンマ関数フィッティングや指数関数による外挿等、任意の公知のアルゴリズムでよい。
【００３３】
続いて、ＣＰＵ５により、動脈入力関数に用いる動脈（以下、流入動脈と記述する）とＰＶＡ効果の補正等の定量化に用いる静脈（以下、流出静脈と記述する）を決定する（ステップ２０５）。流入動脈と流出静脈は、操作者が断層像を見ながら手動で指定してもよく、時間−濃度曲線のピーク値及びピーク時間に基づいて解析対象である器官に対する流入動脈と流出静脈を自動的に決定してもよい。
【００３４】
流入動脈と流出静脈を自動的に決定する場合には、まず、流入動脈と流出静脈の中心画素を算出する。流入動脈と流出静脈の中心画素を自動選択するには、各組織における時間−濃度曲線の最大値と最小値の差（以下、ΔＴＤＣ と記述する）やピーク時間、ピーク値の特性の違いを利用すればよい。例えば、ΔＴＤＣ がある閾値以上の画素の中から、ピーク値が所定の閾値以上で、かつ最もピーク時間の小さい画素を選ぶことにより、流入動脈の中心画素を自動選択することができ、また、ΔＴＤＣ がある閾値以上の画素の中からピーク値が所定の閾値以上で、かつ最もピーク時間の大きい画素を選ぶことにより、流出静脈の中心画素を自動選択することができる。次に、流入動脈と流出静脈の中心画素周辺の画素に対して、ピーク値が閾値以上の画素とそうでない画素を分離し、ピーク値が閾値以上の画素のうち、中心画素を含んだ連結画素を抽出することにより、流入動脈領域と流出静脈領域が決定される。尚、中心画素のみを流入動脈と流出静脈としてもよい。
【００３５】
次に、ＣＰＵ５により、ＰＶＡ効果の補正を行う（ステップ２０６）。尚、ＰＶＡ効果の補正方法については後述する。続いて、ＣＰＵ５により、画素毎に逆フィルタ演算を行い、生体機能情報を算出する（ステップ２０７）。尚、逆フィルタ演算については後述する。
【００３６】
画素毎に解析された生体機能情報をマッピングすることで機能画像が得られ（ステップ２０８）、表示装置８により機能画像が表示される（ステップ２０９）。尚、ステップ２０４は、ステップ２０６とステップ２０７の間で実行してもよい。
【００３７】
次に、不要領域を削除する方法について説明する。
【００３８】
図３は不要領域を削除する方法を示すフローチャートである。同図に示すように、まず、ＣＰＵ５により、図２のステップ２０１で入力した断層像を２値化する（ステップ３０１）。２値化処理は、例えば、断層像において、閾値以上の画素値を有する画素を１に、閾値未満の画素値を有する画素を０にそれぞれ置き換えることで達成される。尚、閾値は、解析対象となる生体組織とルームエアーが分離できるような任意の値でよい。例えば、入力した断層像がＣＴ像であるならば、閾値を−２００程度に設定すればよい。
【００３９】
続いて、ＣＰＵ５により、ステップ３０１で２値化された画像をラベリング処理する（ステップ３０２）。ラベリング処理は、つながっている全ての画素（連結成分）に同じ番号（ラベル）を付けて、異なる連結成分には別の番号を付けることで達成される。尚、隣接する画素が互いに１の場合には、これらの画素はつながっていることになる。ラベルとして付与する値は、例えば５０からの続き番号とする。
【００４０】
次に、ＣＰＵ５により、最大連結成分を探索する（ステップ３０３）。最大連結成分の探索は、ラベリングした画像全体を走査して、ラベル値ごとに画素数を数える処理を行い、最も画素数が多いラベル値を有する連結成分を選び出すことで達成される。その後、ＣＰＵ５により、最大連結成分を抽出する（ステップ３０４）。最大連結成分の抽出は、ステップ３０３で選び出された最も画素数が多いラベル値を有する画素を残し、それ以外の画素を０に置き換えることで達成される。これにより２値化処理の段階では画素値が１に置き換えられる可能性のある寝台や点滴チューブ等の不要領域の画素値が０に置き換えられる。
【００４１】
次に、ＣＰＵ５により、最大連結成分の輪郭を追跡し、最大連結成分の最外周の輪郭線を抽出する（ステップ３０５）。輪郭線の追跡は、ステップ３０４で最大連結成分のみが残された画像の左上隅画素から横に走査し、最初に出会った非ゼロのラベル値を有する画素を開始点として反時計回り方向に輪郭を追跡し、開始点に戻ったら追跡を終了することで達成される。上記輪郭線上の画素にはラベル値とは異なる値、例えば１を代入しておく。
【００４２】
続いて、ＣＰＵ５により、ステップ３０５で求めた輪郭線の内部を１で塗りつぶし、マスク画像を作成する（ステップ３０６）。輪郭線の内部を塗りつぶすには、従来から知られている、閉領域内部の１点を開始点として閉領域内部を塗りつぶす処理であるシードフィルアルゴリズムを適用すればよい。本実施の形態では、閉領域はステップ３０５で求めた輪郭線である。また、閉領域内部の点は、ステップ３０４で求めた最大連結成分のラベル値の付いた画素である。このラベル値の付いた画素を１つ検出し、その点を開始点としてシードフィル処理を行えば、輪郭線内部を塗りつぶすことができる。
【００４３】
次に、ＣＰＵ５により、入力画像（断層像）とステップ３０６で作成したマスク画像とを掛け合わせ、マスク領域外の画素を削除する（ステップ３０７）。具体的には、マスク画像を１画素ずつ調べ、ゼロの画素を検出したら、その画素の座標に対応する原画像の画素を最低輝度に置き換える処理をする。これを全画素について行えば、マスク領域外の画素が入力画像から削除される。
【００４４】
次に、ＰＶＡ効果の補正方法について説明する。
【００４５】
図４は本実施の形態におけるＰＶＡ効果の補正方法を示す図である。例えば、脳の生体機能情報を解析する場合、一般に流入動脈には前大脳動脈や中大脳動脈を、流出静脈には上矢状静脈洞を用いる。Ｌａｐｉｎ らの研究（Ｊ Ｃｏｍｐｕｔ Ａｓｓｉｓｉｔ Ｔｏｍｏｇｒ １９９３；１０８−１１４）によれば、ＰＶＡ効果により動脈を含むボクセルの信号値が低下するため、流入動脈である前大脳動脈や中大脳動脈の時間−濃度曲線は本来の値よりも小さくなる。またＬｅｅ らの研究（ＡＪＮＲ １９９９；２０：６３−７３）によれば、直径１．７３ｍｍ以上の血管はＰＶＡ効果の影響を受けないことから、流出静脈である上矢状静脈洞の時間−濃度曲線はＰＶＡ効果の影響を受けていない。
【００４６】
よって流入動脈と流出静脈の時間濃度曲線は、一般に図４（ａ）のようになり、流入動脈の時間濃度曲線のピーク値Ｐａは、流出静脈のピーク値ＰｖよりもＰＶＡ効果による信号値の低下の分だけ小さな値になる。理想的には血管内の造影剤濃度は血管径の大小に依らず一定であるので、ＰＶＡ効果の影響が無ければ信号値のピーク値は血管径の大小によらず一定である。
【００４７】
従って、流入動脈のピーク値をＰａ、流出静脈のピーク値をＰｖとすると、図４（ｂ）に示すように流入動脈の時間−濃度曲線をＰｖ／Ｐａ倍すれば、流入動脈におけるＰＶＡ効果を補正できる。ＰＶＡ効果を補正する前の時間ｔにおける流入動脈の時間−濃度曲線をＡＩＦｐｖａ（ｔ） 、ＰＶＡ効果を補正した彼の時間ｔにおける流入動脈の時間−濃度曲線をＡＩＦ（ｔ）とすると、両者の間に次式が成り立つように補正すればよい。
【００４８】
【数１】
ＡＩＦ（ｔ）＝ＡＩＦｐｖａ（ｔ） ×Ｐｖ／Ｐａ …（１）
図５は本実施の形態におけるＰＶＡ効果の他の補正方法を示す図である。
【００４９】
過去の臨床画像の中に、これから実施する検査と全く同一の撮影条件、造影条件で撮影した例があるならば、この過去の臨床画像を用いてＰＶＡ効果を補正することも可能である。
【００５０】
例えば、図５（ａ）に示す臨床画像において、断層像内に含まれる複数の血管１〜５に対し、血管径Ｒ１〜Ｒ５とその血管の時間−濃度曲線のピーク値を算出する。これを図表に表すと、図５（ｂ）となり、プロットしてグラフに表すと、図４（ｃ）のようなカーブが得られる。尚、血管径とピーク値の算出方法については後述する。血管径が細いほどＰＶＡ効果の影響を強く受け、血管径がある値以上であればＰＶＡ効果の影響を受けないため、図４（ｃ）のカーブはある径Ｒｔｈまでは徐々に上昇し、Ｒｔｈ以上のＲに対してはピーク値が平衡になる。尚、この平衡状態におけるピーク値を便宜上Ｐｔｈとし、流入動脈における血管径をＲａ、流入動脈における時間−濃度曲線のピーク値をＰａとする。
【００５１】
ＲａがＲｔｈより小さい場合には、流入動脈の時間−濃度曲線をＰｔｈ／Ｐａ倍すれば、流入動脈におけるＰＶＡ効果を補正できる。ＰＶＡ効果を補正する前の時間ｔにおける流入動脈の時間−濃度曲線をＡＩＦｐｖａ（ｔ） 、ＰＶＡ効果を補正した後の時間ｔにおける流入動脈の時間−濃度曲線をＡＩＦ（ｔ）とすると、両者の間に次式が成り立つように補正すればよい。
【００５２】
【数２】
ＡＩＦ（ｔ）＝ＡＩＦｐｖａ（ｔ） ×Ｐｔｈ／Ｐａ …（２）
次に、血管径とピーク値の算出方法について説明する。
【００５３】
図６は本実施の形態における血管径とピーク値の算出手順を示すフローチャートである。まず、ＣＰＵ５により、各時間−濃度曲線のピーク時間とピーク値を算出する（ステップ６０１）。次に、ＣＰＵ５により、血管の中心画素を算出する（ステップ６０２）。中心画素を算出するには、ピーク値の特性を利用すればよく、ピーク値が最大の画素を選ぶことにより、血管の中心画素を自動選択することができる。続いて、ＣＰＵ５により、血管の中心画素周辺の画素に対して、ピーク値が閾値以上の画素とそうでない画素を分離する（ステップ６０３）。その後、ステップ６０３で求めたピーク値が閾値以上の画素のうち、中心画素を含んだ連結画素を抽出することにより、血管領域が決定される（ステップ６０４）。
【００５４】
次に、ＣＰＵ５により、血管領域内の画素数を算出し（ステップ６０５）、その算出した画素数から血管領域の面積を算出する（ステップ６０６）。一画素の面積は装置毎に一定値であるので、一画素あたりの面積にステップ６０５で求められた画素数を掛けることで血管領域の面積Ｓが得られる。
【００５５】
続いて、ＣＰＵ５により、血管径Ｒを算出する（ステップ６０７）。断層像上での血管の形状は円形で近似して差し支えないので、血管径Ｒと面積Ｓには、次式、
【００５６】
【数３】
Ｒ＝√（Ｓ／π） …（３）
が成り立ち、（３）式に従って血管径を算出すればよい。
【００５７】
次に、ＣＰＵ５により血管の時間−濃度曲線を算出する（ステップ６０８）。ステップ２０３において既に各画素における時間−濃度曲線が算出されているので、血管領域の全ての画素における時間−濃度曲線の平均値を算出することで、血管の時間−濃度曲線が得られる。また、平均値をとるかわりに、ステップ６０２で算出されている血管中心画素からの距離に応じて、任意に重み付けをしてもよい。
【００５８】
次に、ＣＰＵ５により、血管の時間−濃度曲線のピーク値を算出する（ステップ６０９）。尚、ステップ６０５からステップ６０７までの処理と、ステップ６０８からステップ６０９までの処理は、処理の順序が逆になってもよい。本実施例では、ピーク値の特性を利用して血管の抽出を行ったが、血管の抽出法はこれに限定されるものではなく、血管が正確に抽出できる方法であるならば任意の方法でよい。また、断層像を見ながら手動で血管領域を指定することにより、血管を抽出してもよい。
【００５９】
次に、逆フィルタ法の概念について説明する。
【００６０】
図７は本実施の形態における逆フィルタ法を説明するための概念図である。伝達関数ｆの対象組織において、入力関数をＣｉｎ（ｔ）、出力関数をＣｏｕｔ（ｔ） とおくと、Ｃｏｕｔ（ｔ） はＣｉｎ（ｔ）とｆを畳み込み積分したものになり、次式、
【００６１】
【数４】
のように表される。
【００６２】
造影剤を動脈から体内に注入する場合、Ｃｉｎ（ｔ）は流入動脈における時間−濃度曲線（動脈入力関数）ＡＩＦ（ｔ）に、Ｃｏｕｔ（ｔ） は組織における時間−濃度曲線Ｃｒ（ｔ） に置き換えられ、次式、
【００６３】
【数５】
のように表される。
【００６４】
（５）式を周波数空間に変換すると、次式、
【００６５】
【数６】
Ｆ｛Ｃｒ（ｔ）｝＝Ｆ｛ＡＩＦ（ｔ）｝・Ｆ｛ｆ（ｔ）｝…（６）
のようになる。
【００６６】
ここで、Ｆ｛ＡＩＦ（ｔ）｝はＡＩＦ（ｔ）のフーリエ変換を、Ｆ｛Ｃｒ（ｔ） ｝はＣｒ（ｔ） のフーリエ変換を、Ｆ｛ｆ（ｔ）｝はｆ（ｔ）のフーリエ変換をそれぞれ表す。
【００６７】
（５）式と（６）式より、流入動脈における時間−濃度曲線と組織における時間−濃度曲線を逆畳み込み演算すれば伝達関数ｆを求めることができるが、逆畳み込み演算を行う代わりに、流入動脈における時間−濃度曲線の逆フィルタを組織の時間−濃度曲線に掛けても等価な演算結果が得られる。
【００６８】
ＡＩＦ（ｔ）の逆フィルタをＡＩＦ（ｔ）−１とおくと、（６）式より、ｆ（ｔ）は次式、
【００６９】
【数７】
ｆ（ｔ）＝Ｆ−１［Ｆ｛ＡＩＦ−１（ｔ） ｝・Ｆ｛Ｃｒ（ｔ） ｝］ …（７）
のように表される。
【００７０】
ここでＦ−１は、逆フーリエ変換を表す。
【００７１】
図８は本実施の形態における逆フィルタ法を説明するための概念図である。
【００７２】
血流動態解析を行う場合、被検者への侵襲性を抑えるため、たいていの場合、造影剤は静脈から体内へ注入する。この場合、造影剤は体内を循環した後、解析対象組織に流入する。入力関数をＣｉｎ（ｔ）、循環経路における伝達関数をｆ１ 、対象組織における伝達関数をｆ２ 、対象組織からの出力関数をＣｏｕｔ（ｔ） とおくと、Ｃｏｕｔ（ｔ） は、Ｃｉｎ（ｔ）とｆ１ を畳み込み積分したものを、更にｆ２ と畳み込み積分したものになるので、次式、
【００７３】
【数８】
のようになる。
【００７４】
Ｃｉｎ（ｔ）とｆ１ のたたみこみ積分は流入動脈における時間−濃度曲線（動脈入力関数）ＡＩＦ（ｔ）に、Ｃｏｕｔ（ｔ） は組織における時間−濃度曲線Ｃｒ（ｔ） に置き換えられ、次式、
【００７５】
【数９】
のように表される。
【００７６】
（９）式は、前述した（５）式と等価である。よって伝達関数ｆ２ は、次式、
【００７７】
【数１０】
ｆ２（ｔ） ＝Ｆ−１［Ｆ｛ＡＩＦ−１（ｔ） ｝・Ｆ｛Ｃｒ（ｔ） ｝］ …（１０）
のように表される。
【００７８】
次に、逆フィルタ演算について説明する。
【００７９】
図９は本実施の形態における逆フィルタ演算の処理手順を示すフローチャートである。同図に示すように、ＣＰＵ５により、流入動脈の時間−濃度曲線をフーリエ変換する（ステップ９０１）。次に、ＣＰＵ５により、流入動脈の時間−濃度曲線のフーリエ変換から逆フィルタを算出する（ステップ９０２）。流入動脈の時間−濃度曲線のフーリエ変換の実部をＦ｛ＡＩＦＲ｝、虚部をＦ｛ＡＩＦＩ｝とすると、逆フィルタは、次式、
【００８０】
【数１１】
で表される。
【００８１】
次に、逆フィルタの高周波数成分を制御する（ステップ９０３）。ウィナーフィルタやバターワースフィルタなどの高周波数領域制御フィルタによって逆フィルタをフィルタリングすることにより、逆フィルタの高周波数成分が制御される。
【００８２】
続いて、ＣＰＵ５により、各組織の時間−濃度曲線をフーリエ変換する（ステップ９０４）。このフーリエ変換に対して、ステップ９０２で求めた逆フィルタでフィルタリングする（ステップ９０５）。このフィルタリングされた組織の時間−濃度曲線のフーリエ変換を逆フーリエ変換する伝達関数が求められる（ステップ９０６）。
【００８３】
次に、ＣＰＵ５により、ステップ９０６で求められた伝達関数から、生体機能情報を算出する（ステップ９０７）。生体機能情報の算出方法については後述する。尚、ステップ９０４からステップ９０７までは、解析する組織の断層像上での画素数分だけ繰り返す。また、ステップ９０１からステップ９０３までと、ステップ９０４は、処理の順序は逆になってもよい。
【００８４】
次に、伝達関数から生体機能情報を算出する方法について説明する。
【００８５】
図１０は生体機能情報の算出方法の説明図である。伝達関数をｆ（ｔ）とすると、ｆ（ｔ）は、例えば図１０のような形状になり、血流量（ｂｌｏｏｄｆｌｏｗ：以下、ＢＦと記述する）は、次式、
【００８６】
【数１２】
ＢＦ＝ｆｍａｘ …（１２）
のようにｆ（ｔ）の最大値から求められ、血液量（ｂｌｏｏｄｖｏｌｕｍｅ：以下、ＢＶと記述する）は、次式、
【００８７】
【数１３】
ＢＶ＝∫ｆ（ｔ）ｄｔ …（１３）
のようにｆ（ｔ）の曲線下面積から求められ、平均通過時間（ｍｅａｎｔｒａｎｓｉｔ ｔｉｍｅ： 以下、ＭＴＴ と記述する）は、次式、
【００８８】
【数１４】
ＭＴＴ ＝∫ｆ（ｔ）ｄｔ／ｆｍａｘ …（１４）
のようにｆ（ｔ）の幅から求められる。
【００８９】
ここで、ｆ ｍａｘ は、ｆ（ｔ）の最大値を表す。
【００９０】
【発明の効果】
以上説明したように本発明によれば、コンピュータ断層診断装置の提供するダイナミック撮影した断層像から生体機能情報解析を行う場合、ガンマ関数解析法や最大勾配法よりも低い造影レートで、ｄｅｃｏｎｖｏｌｕｔｉｏｎ 法よりも短い演算時間で、生体機能情報解析を行えるという効果がある。また、ＰＶＡ効果による定量性の悪化を防止できるという効果もある。更に、操作者への負担を増加させることなく不要領域の影響による無用な演算時間の増加や解析精度の悪化を防止する事ができるという効果もある。
【図面の簡単な説明】
【図１】本発明に係る血流動態解析装置のハードウェア構成例を示すブロック図
【図２】本発明に係る血流動態解析装置の断層像の入力から機能画像の表示までを示すフローチャート
【図３】不要領域を削除する方法を示すフローチャート
【図４】本実施の形態におけるＰＶＡ効果の補正方法を示す図
【図５】本実施の形態におけるＰＶＡ効果の他の補正方法を示す図
【図６】本実施の形態における血管径とピーク値の算出手順を示すフローチャート
【図７】本実施の形態における逆フィルタ法を説明するための概念図
【図８】本実施の形態における逆フィルタ法を説明するための概念図
【図９】本実施の形態における逆フィルタ演算の処理手順を示すフローチャート
【図１０】生体機能情報の算出方法の説明図
【符号の説明】
１…断層像入力装置、２…演算装置、３…インターフェイス（Ｉ／Ｆ）、４…モリ、５…中央処理装置（ＣＰＵ）、６…記録媒体、７…外部入力装置、８…表示装置、９…共通バス
【発明の属する技術分野】
本発明は血流動態解析装置に係り、特にＸ線や電磁波を用いて被検体の時間変化に依存する生体機能情報を得ることが可能なコンピュータ断層診断装置の提供する断層像から生体機能情報解析等の画像解析を行うことに対して有用な技術に関する。
【０００２】
【従来の技術】
従来、ダイナミック撮影を行うコンピュータ断層診断装置として、陽電子放射断層撮影装置（Ｐｏｓｉｔｒｏｎ Ｅｍｉｓｓｉｏｎ Ｔｏｍｏｇｒａｐｈｙ：ＰＥＴ）や単光子放射ＣＴ（ＳｉｎｇｌｅＰｈｏｔｏｎ Ｅｍｉｓｓｉｏｎ ＣＴ：ＳＰＥＣＴ）などの核医学装置がある。核医学装置による血流動態計測は、放射性核種を被検体内に注入し、核種からの放射線をシンチレーションカメラで計測すると、横断面の放射性核種分布像を得られ、この放射性核種分布像を解析することで、器官の生体機能情報を解析していた。
【０００３】
他のダイナミック撮影を行うコンピュータ断層診断装置として、Ｘ線ＣＴ装置がある。単純Ｘ線ＣＴ像では病変の超早期相の診断が困難である場合には、Ｘ線ＣＴ装置によるダイナミック撮影を行い、この撮影時にコントラストを付けるために主にヨード系の造影剤を用いていた。
【０００４】
造影剤を注入してダイナミックスキャンを行うと、造影剤濃度とＣＴ値は比例関係にあることから、撮影断面の時間変化の情報が得られる。撮影断面の時間変化から各組織のＣＴ値の時間変化である時間−濃度曲線が得られ、この時間−濃度曲線を解析することで器官の生体機能情報を解析していた（例えば非特許文献１参照。）。
【０００５】
更に他のダイナミック撮影を行うコンピュータ断層診断装置として、ＭＲ装置がある。高磁場ＭＲ装置では非造影の診断方法もあるが、高磁場ＭＲ装置でコントラストの良い画像を得たい場合や低磁場ＭＲ装置においては、主にガドリニウム系の造影剤を用いていた。ＭＲ装置においても撮影断面の時間変化から各組織のＭＲ値の相対的信号強度変化率である時間−濃度曲線が得られる。この時間−濃度曲線を解析することで器官の時間変化に依存する生体機能情報を解析していた（例えば非特許文献２参照。）。
【０００６】
一方、血流動態を解析するアルゴリズムとして、ガンマ関数解析法がある。ガンマ関数解析法では、時間−濃度曲線をガンマ関数でフィッティングし、フィッティングカーブのピーク値や曲線下面積から血流情報を算出していた（例えば非特許文献３参照。）。
【０００７】
他の血流動態を解析するアルゴリズムとして、最大勾配法（ｍａｘｉｍｕｍ ｓｌｏｐｅ 法）がある。最大勾配法では、各組織における時間−濃度曲線の傾きの最大値を動脈入力関数におけるＣＴ値上昇の最大値で割ることで血流量を算出していた（例えば非特許文献４参照。）。
【０００８】
更に他の血流動態を解析するアルゴリズムとして、ｄｅｃｏｎｖｏｌｕｔｉｏｎ（デコンボリューション）法がある。ｄｅｃｏｎｖｏｌｕｔｉｏｎ 法では、動脈入力関数と組織残留関数をｄｅｃｏｎｖｏｌｕｔｉｏｎ 演算することでインパルス残留関数を求め、インパルス残留関数のピーク値や曲線下面積から血流情報を算出していた（例えば非特許文献５参照。）。
【０００９】
【非特許文献１】
Ｊ．Ｐ．Ｍｕｉｚｅｌａａｒ ｅｔｃ．：Ａ Ｎｅｗ Ｍｅｔｈｏｄ ｆｏｒ Ｑｕａｎｔｉｔａｔｉｖｅ Ｒｅｇｉｏｎａｌ Ｃｅｒｅｂｒａｌ ＢｌｏｏｄＶｏｌｕｍｅ Ｍｅａｓｕｒｅｍｅｎｔｓ Ｕｓｉｎｇ Ｃｏｍｐｕｔｅｄ Ｔｏｍｏｇｒａｐｈｙ．：Ｓｔｒｏｋｅ１９９７；２８：１９９８−２００５
【００１０】
【非特許文献２】
井川房雄 他：拡散強調画像と潅流画像による急性期脳梗塞の診断の意義：日磁医誌 第２０巻５号２０３−２０９頁（２０００年）
【００１１】
【非特許文献３】
Ｌ．Ａｘｅｌ：Ｃｅｒｅｂｒａｌ Ｂｌｏｏｄ Ｆｌｏｗ Ｄｅｔｅｒｍｉｎａｔｉｏｎ ｂｙ Ｒａｐｉｄ−Ｓｅｑｕｅｎｃｅ Ｃｏｍｐｕｔｅｄ Ｔｏｍｏｇｒａｐｈｙ：Ｒａｄｉｏｌｏｇｙ １９８０；１３７：６７９−６８６
【００１２】
【非特許文献４】
Ｍ．Ｋｏｅｎｉｇ ｅｔｃ．：Ｐｅｒｆｕｓｉｏｎ ＣＴ ｏｆ ｔｈｅ Ｂｒａｉｎ：Ｒａｄｉｏｌｏｇｙ
１９９８；２０９：８５−９３
【００１３】
【非特許文献５】
Ｌ．Ｏｓｔｅｒｇａａｒｄ ｅｔｃ．：Ｈｉｇｈ Ｒｅｓｏｌｕｔｉｏｎ Ｍｅａｓｕｒｅｍｅｎｔ ｏｆ Ｃｅｒｅｂｒａｌ Ｂｌｏｏｄ Ｆｌｏｗ ｕｓｉｎｇ Ｉｎｔｒａｖａｓｃｕｌａｒ Ｔｒａｃｅｒ Ｂｏｌｕｓ Ｐａｓｓａｇｅｓ：１９９６；ＮＲＭ ３６：７１５−７２５
【００１４】
【発明が解決しようとする課題】
本発明者は、上記従来技術を検討した結果、以下の問題点を見い出した。ガンマ関数解析法は、８〜１０ｍｌ／秒程度の高い造影レートが必要であるため侵襲性が高く、定量性も悪いという問題があり、同様に最大勾配法は、８ｍｌ／秒程度の高い造影レートが必要であるため侵襲性が高いという問題がある。
【００１５】
一方、ｄｅｃｏｎｖｏｌｕｔｉｏｎ 法は、３〜５ｍｌ／秒程度の低い造影レートで検査可能であるため侵襲性が低く、定量性も優れているが、動脈入力関数と組織残留関数のｄｅｃｏｎｖｏｌｕｔｉｏｎ 演算は非常に演算時間がかかるため、解析結果が出るまでに長時間必要であるという問題がある。特に、急性期脳梗塞のように迅速性を求められる症例においては、血流異常の有無の診断は出来る限り早急に行わなければならないため、演算時間の長さが問題になる。
【００１６】
また、ＣＴ装置を用いて血流動態を検査する場合、パーシャルボリュームアベレージング（以下、ＰＶＡと記述する）効果により動脈を含むボクセルのＣＴ値は本来の値よりも低下する。よって動脈入力関数が本来の値よりも小さくなるため、解析結果の定量性が悪化するという問題がある。
【００１７】
更に、断層像中には、ルームエアーや寝台、骨等の生体機能情報の解析には不要な部分も存在する。断層像中に不要な部分が存在すると、不要領域の信号値の影響で解析精度が悪化するという問題がある。また、不要領域を削除しない場合、不要領域の分だけ無用に演算時間が長くなるという問題もある。更に、不要領域を手動で除去しようとすると、多くの作業時間が消費され、操作者への負担が増大するという問題もある。
【００１８】
本発明の目的は、コンピュータ断層診断装置の提供する断層像から解析対象とする器官の生体機能情報を得る場合において、低い造影レートで定量性のよい生体機能情報を得る事ができ、特に出来る限り短い演算時間で生体機能情報を得る事が可能な血流動態解析装置を提供することにある。
【００１９】
本発明の他の目的は、コンピュータ断層診断装置の提供する断層像から解析対象とする器官の生体機能情報を得る場合において、操作者への負担を増加させることなく解析に不要な領域の信号値の影響による解析精度の悪化を防止する事が可能な血流動態解析装置を提供することにある。
【００２０】
本発明の更に他の目的は、コンピュータ断層診断装置の提供する断層像から解析対象とする器官の生体機能情報を得る場合において、ＰＶＡ効果の影響による定量性の悪化を防止する事が可能な血流動態解析装置を提供することにある。
【００２１】
【課題を解決するための手段】
上記問題を解決するために請求項１に係る血流動態解析装置は、コンピュータ断層像撮影手段でダイナミック撮影した断層像の入力手段と、前記入力手段から入力した断層像に基づいて各画素ごとに画素値の時間変化の情報を表す時間−濃度曲線を求める演算手段であって、流入動脈の時間−濃度曲線及び各組織における時間−濃度曲線を求める演算手段と、前記求めた流入動脈の時間−濃度曲線の逆フィルタを算出する手段と、前記算出した逆フィルタを前記求めた各組織における時間−濃度曲線に掛け合わせることで各組織の伝達関数を求める手段と、各組織の伝達関数を使用して生体機能情報を算出する手段と、を備えたことを特徴としている。
【００２２】
即ち、流入動脈の時間−濃度曲線の逆フィルタと各組織における時間−濃度曲線との周波数空間上での掛け算により各組織の伝達関数を求めることができ、従来のｄｅｃｏｎｖｏｌｕｔｉｏｎ 法よりも演算時間の短縮化を図ることができる。
【００２３】
請求項２に示すように前記演算手段は、断層像内の最大連結画素領域を抽出する手段と、前記抽出した最大連結画素領域を示すマスク画像を作成する手段と、前記マスク画像を使用して前記入力手段から入力する断層像のうちの当該マスク画像の領域以外の画素を削除する手段とを有し、ルームエアーや寝台、骨等の生体機能情報解析には不要な領域を自動的に除去することを特徴としている。
【００２４】
請求項３に示すように前記演算手段は、流出静脈の時間−濃度曲線を求めるとともに、前記求めた流入動脈の時間−濃度曲線及び流出静脈の時間−濃度曲線のピーク値をそれぞれ求める手段と、前記流入動脈の時間−濃度曲線のピーク値が前記流出静脈の時間−濃度曲線のピーク値が一致するように前記流入動脈の時間−濃度曲線を補正する補正手段とを有し、流入動脈の時間−濃度曲線におけるパーシャルボリュームアベレージング効果を補正することを特徴としている。
【００２５】
【発明の実施の形態】
以下、添付図面に従って本発明に係る血流動態解析装置の好ましい実施の形態について詳説する。
【００２６】
図１は本発明に係る血流動態解析装置のハードウェア構成例を示すブロック図である。同図に示すように、この血流動態解析装置は、主としてコンピュータ断層像撮影手段で撮影した断層像を入力する断層像入力装置１と、画像解析などの諸演算を行う演算装置２とから構成されている。
【００２７】
断層像入力装置１と演算装置２は必ずしも別個の装置である必要はなく、一体化していてもよい。また、断層像入力装置１と演算装置２は、必ずしもコンピュータ断層像撮影手段と別個の装置である必要はなく、一体化していてもよい。
【００２８】
演算装置２は、断層像入力装置１から断層像の取り込みを行うインターフェイス（Ｉ／Ｆ）３と、断層像や演算結果を一時的に格納するメモリ４と、諸演算を行う中央処理装置（ＣＰＵ）５と、演算結果や機能画像を記録するハードディスクなどの記録媒体６と、マウスやキーボードなどの外部入力装置７と、処理結果を表示する表示装置８とから構成され、これらは共通バス９によって接続されている。
【００２９】
図２は本発明に係る血流動態解析装置の断層像の入力から機能画像の表示までを示すフローチャートである。
【００３０】
まず、断層像入力装置１から血流動態解析を行う被検体の断層像を選択し、選択された断層像を演算装置２に入力する（ステップ２０１）。入力された断層像はメモリ４に一時格納されるか、または記録媒体６に保存される。断層像入力装置１と演算装置２とが一体化している場合には、外部入力装置７により血流動態解析を行う被検体の断層像を選択し、選択された断層像を記録媒体６から読み出し、メモリ４に一時格納する。
【００３１】
続いて、ＣＰＵ５により、断層像に写り込んでいるルームエアーや寝台、骨等、生体機能情報の解析には不要な領域を取り除く（ステップ２０２）。不要領域除去の方法については後述する。
【００３２】
次に、ＣＰＵ５により、断層像から画素ごとに画素値の時間変化の情報を表す時間−濃度曲線を取り出す（ステップ２０３）。その後、ＣＰＵ５により、時間−濃度曲線から第一循環成分を取り出す（ステップ２０４）。第一循環成分を取り出すアルゴリズムは、ガンマ関数フィッティングや指数関数による外挿等、任意の公知のアルゴリズムでよい。
【００３３】
続いて、ＣＰＵ５により、動脈入力関数に用いる動脈（以下、流入動脈と記述する）とＰＶＡ効果の補正等の定量化に用いる静脈（以下、流出静脈と記述する）を決定する（ステップ２０５）。流入動脈と流出静脈は、操作者が断層像を見ながら手動で指定してもよく、時間−濃度曲線のピーク値及びピーク時間に基づいて解析対象である器官に対する流入動脈と流出静脈を自動的に決定してもよい。
【００３４】
流入動脈と流出静脈を自動的に決定する場合には、まず、流入動脈と流出静脈の中心画素を算出する。流入動脈と流出静脈の中心画素を自動選択するには、各組織における時間−濃度曲線の最大値と最小値の差（以下、ΔＴＤＣ と記述する）やピーク時間、ピーク値の特性の違いを利用すればよい。例えば、ΔＴＤＣ がある閾値以上の画素の中から、ピーク値が所定の閾値以上で、かつ最もピーク時間の小さい画素を選ぶことにより、流入動脈の中心画素を自動選択することができ、また、ΔＴＤＣ がある閾値以上の画素の中からピーク値が所定の閾値以上で、かつ最もピーク時間の大きい画素を選ぶことにより、流出静脈の中心画素を自動選択することができる。次に、流入動脈と流出静脈の中心画素周辺の画素に対して、ピーク値が閾値以上の画素とそうでない画素を分離し、ピーク値が閾値以上の画素のうち、中心画素を含んだ連結画素を抽出することにより、流入動脈領域と流出静脈領域が決定される。尚、中心画素のみを流入動脈と流出静脈としてもよい。
【００３５】
次に、ＣＰＵ５により、ＰＶＡ効果の補正を行う（ステップ２０６）。尚、ＰＶＡ効果の補正方法については後述する。続いて、ＣＰＵ５により、画素毎に逆フィルタ演算を行い、生体機能情報を算出する（ステップ２０７）。尚、逆フィルタ演算については後述する。
【００３６】
画素毎に解析された生体機能情報をマッピングすることで機能画像が得られ（ステップ２０８）、表示装置８により機能画像が表示される（ステップ２０９）。尚、ステップ２０４は、ステップ２０６とステップ２０７の間で実行してもよい。
【００３７】
次に、不要領域を削除する方法について説明する。
【００３８】
図３は不要領域を削除する方法を示すフローチャートである。同図に示すように、まず、ＣＰＵ５により、図２のステップ２０１で入力した断層像を２値化する（ステップ３０１）。２値化処理は、例えば、断層像において、閾値以上の画素値を有する画素を１に、閾値未満の画素値を有する画素を０にそれぞれ置き換えることで達成される。尚、閾値は、解析対象となる生体組織とルームエアーが分離できるような任意の値でよい。例えば、入力した断層像がＣＴ像であるならば、閾値を−２００程度に設定すればよい。
【００３９】
続いて、ＣＰＵ５により、ステップ３０１で２値化された画像をラベリング処理する（ステップ３０２）。ラベリング処理は、つながっている全ての画素（連結成分）に同じ番号（ラベル）を付けて、異なる連結成分には別の番号を付けることで達成される。尚、隣接する画素が互いに１の場合には、これらの画素はつながっていることになる。ラベルとして付与する値は、例えば５０からの続き番号とする。
【００４０】
次に、ＣＰＵ５により、最大連結成分を探索する（ステップ３０３）。最大連結成分の探索は、ラベリングした画像全体を走査して、ラベル値ごとに画素数を数える処理を行い、最も画素数が多いラベル値を有する連結成分を選び出すことで達成される。その後、ＣＰＵ５により、最大連結成分を抽出する（ステップ３０４）。最大連結成分の抽出は、ステップ３０３で選び出された最も画素数が多いラベル値を有する画素を残し、それ以外の画素を０に置き換えることで達成される。これにより２値化処理の段階では画素値が１に置き換えられる可能性のある寝台や点滴チューブ等の不要領域の画素値が０に置き換えられる。
【００４１】
次に、ＣＰＵ５により、最大連結成分の輪郭を追跡し、最大連結成分の最外周の輪郭線を抽出する（ステップ３０５）。輪郭線の追跡は、ステップ３０４で最大連結成分のみが残された画像の左上隅画素から横に走査し、最初に出会った非ゼロのラベル値を有する画素を開始点として反時計回り方向に輪郭を追跡し、開始点に戻ったら追跡を終了することで達成される。上記輪郭線上の画素にはラベル値とは異なる値、例えば１を代入しておく。
【００４２】
続いて、ＣＰＵ５により、ステップ３０５で求めた輪郭線の内部を１で塗りつぶし、マスク画像を作成する（ステップ３０６）。輪郭線の内部を塗りつぶすには、従来から知られている、閉領域内部の１点を開始点として閉領域内部を塗りつぶす処理であるシードフィルアルゴリズムを適用すればよい。本実施の形態では、閉領域はステップ３０５で求めた輪郭線である。また、閉領域内部の点は、ステップ３０４で求めた最大連結成分のラベル値の付いた画素である。このラベル値の付いた画素を１つ検出し、その点を開始点としてシードフィル処理を行えば、輪郭線内部を塗りつぶすことができる。
【００４３】
次に、ＣＰＵ５により、入力画像（断層像）とステップ３０６で作成したマスク画像とを掛け合わせ、マスク領域外の画素を削除する（ステップ３０７）。具体的には、マスク画像を１画素ずつ調べ、ゼロの画素を検出したら、その画素の座標に対応する原画像の画素を最低輝度に置き換える処理をする。これを全画素について行えば、マスク領域外の画素が入力画像から削除される。
【００４４】
次に、ＰＶＡ効果の補正方法について説明する。
【００４５】
図４は本実施の形態におけるＰＶＡ効果の補正方法を示す図である。例えば、脳の生体機能情報を解析する場合、一般に流入動脈には前大脳動脈や中大脳動脈を、流出静脈には上矢状静脈洞を用いる。Ｌａｐｉｎ らの研究（Ｊ Ｃｏｍｐｕｔ Ａｓｓｉｓｉｔ Ｔｏｍｏｇｒ １９９３；１０８−１１４）によれば、ＰＶＡ効果により動脈を含むボクセルの信号値が低下するため、流入動脈である前大脳動脈や中大脳動脈の時間−濃度曲線は本来の値よりも小さくなる。またＬｅｅ らの研究（ＡＪＮＲ １９９９；２０：６３−７３）によれば、直径１．７３ｍｍ以上の血管はＰＶＡ効果の影響を受けないことから、流出静脈である上矢状静脈洞の時間−濃度曲線はＰＶＡ効果の影響を受けていない。
【００４６】
よって流入動脈と流出静脈の時間濃度曲線は、一般に図４（ａ）のようになり、流入動脈の時間濃度曲線のピーク値Ｐａは、流出静脈のピーク値ＰｖよりもＰＶＡ効果による信号値の低下の分だけ小さな値になる。理想的には血管内の造影剤濃度は血管径の大小に依らず一定であるので、ＰＶＡ効果の影響が無ければ信号値のピーク値は血管径の大小によらず一定である。
【００４７】
従って、流入動脈のピーク値をＰａ、流出静脈のピーク値をＰｖとすると、図４（ｂ）に示すように流入動脈の時間−濃度曲線をＰｖ／Ｐａ倍すれば、流入動脈におけるＰＶＡ効果を補正できる。ＰＶＡ効果を補正する前の時間ｔにおける流入動脈の時間−濃度曲線をＡＩＦｐｖａ（ｔ） 、ＰＶＡ効果を補正した彼の時間ｔにおける流入動脈の時間−濃度曲線をＡＩＦ（ｔ）とすると、両者の間に次式が成り立つように補正すればよい。
【００４８】
【数１】
ＡＩＦ（ｔ）＝ＡＩＦｐｖａ（ｔ） ×Ｐｖ／Ｐａ …（１）
図５は本実施の形態におけるＰＶＡ効果の他の補正方法を示す図である。
【００４９】
過去の臨床画像の中に、これから実施する検査と全く同一の撮影条件、造影条件で撮影した例があるならば、この過去の臨床画像を用いてＰＶＡ効果を補正することも可能である。
【００５０】
例えば、図５（ａ）に示す臨床画像において、断層像内に含まれる複数の血管１〜５に対し、血管径Ｒ１〜Ｒ５とその血管の時間−濃度曲線のピーク値を算出する。これを図表に表すと、図５（ｂ）となり、プロットしてグラフに表すと、図４（ｃ）のようなカーブが得られる。尚、血管径とピーク値の算出方法については後述する。血管径が細いほどＰＶＡ効果の影響を強く受け、血管径がある値以上であればＰＶＡ効果の影響を受けないため、図４（ｃ）のカーブはある径Ｒｔｈまでは徐々に上昇し、Ｒｔｈ以上のＲに対してはピーク値が平衡になる。尚、この平衡状態におけるピーク値を便宜上Ｐｔｈとし、流入動脈における血管径をＲａ、流入動脈における時間−濃度曲線のピーク値をＰａとする。
【００５１】
ＲａがＲｔｈより小さい場合には、流入動脈の時間−濃度曲線をＰｔｈ／Ｐａ倍すれば、流入動脈におけるＰＶＡ効果を補正できる。ＰＶＡ効果を補正する前の時間ｔにおける流入動脈の時間−濃度曲線をＡＩＦｐｖａ（ｔ） 、ＰＶＡ効果を補正した後の時間ｔにおける流入動脈の時間−濃度曲線をＡＩＦ（ｔ）とすると、両者の間に次式が成り立つように補正すればよい。
【００５２】
【数２】
ＡＩＦ（ｔ）＝ＡＩＦｐｖａ（ｔ） ×Ｐｔｈ／Ｐａ …（２）
次に、血管径とピーク値の算出方法について説明する。
【００５３】
図６は本実施の形態における血管径とピーク値の算出手順を示すフローチャートである。まず、ＣＰＵ５により、各時間−濃度曲線のピーク時間とピーク値を算出する（ステップ６０１）。次に、ＣＰＵ５により、血管の中心画素を算出する（ステップ６０２）。中心画素を算出するには、ピーク値の特性を利用すればよく、ピーク値が最大の画素を選ぶことにより、血管の中心画素を自動選択することができる。続いて、ＣＰＵ５により、血管の中心画素周辺の画素に対して、ピーク値が閾値以上の画素とそうでない画素を分離する（ステップ６０３）。その後、ステップ６０３で求めたピーク値が閾値以上の画素のうち、中心画素を含んだ連結画素を抽出することにより、血管領域が決定される（ステップ６０４）。
【００５４】
次に、ＣＰＵ５により、血管領域内の画素数を算出し（ステップ６０５）、その算出した画素数から血管領域の面積を算出する（ステップ６０６）。一画素の面積は装置毎に一定値であるので、一画素あたりの面積にステップ６０５で求められた画素数を掛けることで血管領域の面積Ｓが得られる。
【００５５】
続いて、ＣＰＵ５により、血管径Ｒを算出する（ステップ６０７）。断層像上での血管の形状は円形で近似して差し支えないので、血管径Ｒと面積Ｓには、次式、
【００５６】
【数３】
Ｒ＝√（Ｓ／π） …（３）
が成り立ち、（３）式に従って血管径を算出すればよい。
【００５７】
次に、ＣＰＵ５により血管の時間−濃度曲線を算出する（ステップ６０８）。ステップ２０３において既に各画素における時間−濃度曲線が算出されているので、血管領域の全ての画素における時間−濃度曲線の平均値を算出することで、血管の時間−濃度曲線が得られる。また、平均値をとるかわりに、ステップ６０２で算出されている血管中心画素からの距離に応じて、任意に重み付けをしてもよい。
【００５８】
次に、ＣＰＵ５により、血管の時間−濃度曲線のピーク値を算出する（ステップ６０９）。尚、ステップ６０５からステップ６０７までの処理と、ステップ６０８からステップ６０９までの処理は、処理の順序が逆になってもよい。本実施例では、ピーク値の特性を利用して血管の抽出を行ったが、血管の抽出法はこれに限定されるものではなく、血管が正確に抽出できる方法であるならば任意の方法でよい。また、断層像を見ながら手動で血管領域を指定することにより、血管を抽出してもよい。
【００５９】
次に、逆フィルタ法の概念について説明する。
【００６０】
図７は本実施の形態における逆フィルタ法を説明するための概念図である。伝達関数ｆの対象組織において、入力関数をＣｉｎ（ｔ）、出力関数をＣｏｕｔ（ｔ） とおくと、Ｃｏｕｔ（ｔ） はＣｉｎ（ｔ）とｆを畳み込み積分したものになり、次式、
【００６１】
【数４】
【００６２】
造影剤を動脈から体内に注入する場合、Ｃｉｎ（ｔ）は流入動脈における時間−濃度曲線（動脈入力関数）ＡＩＦ（ｔ）に、Ｃｏｕｔ（ｔ） は組織における時間−濃度曲線Ｃｒ（ｔ） に置き換えられ、次式、
【００６３】
【数５】
【００６４】
（５）式を周波数空間に変換すると、次式、
【００６５】
【数６】
Ｆ｛Ｃｒ（ｔ）｝＝Ｆ｛ＡＩＦ（ｔ）｝・Ｆ｛ｆ（ｔ）｝…（６）
のようになる。
【００６６】
ここで、Ｆ｛ＡＩＦ（ｔ）｝はＡＩＦ（ｔ）のフーリエ変換を、Ｆ｛Ｃｒ（ｔ） ｝はＣｒ（ｔ） のフーリエ変換を、Ｆ｛ｆ（ｔ）｝はｆ（ｔ）のフーリエ変換をそれぞれ表す。
【００６７】
（５）式と（６）式より、流入動脈における時間−濃度曲線と組織における時間−濃度曲線を逆畳み込み演算すれば伝達関数ｆを求めることができるが、逆畳み込み演算を行う代わりに、流入動脈における時間−濃度曲線の逆フィルタを組織の時間−濃度曲線に掛けても等価な演算結果が得られる。
【００６８】
ＡＩＦ（ｔ）の逆フィルタをＡＩＦ（ｔ）−１とおくと、（６）式より、ｆ（ｔ）は次式、
【００６９】
【数７】
ｆ（ｔ）＝Ｆ−１［Ｆ｛ＡＩＦ−１（ｔ） ｝・Ｆ｛Ｃｒ（ｔ） ｝］ …（７）
のように表される。
【００７０】
ここでＦ−１は、逆フーリエ変換を表す。
【００７１】
図８は本実施の形態における逆フィルタ法を説明するための概念図である。
【００７２】
血流動態解析を行う場合、被検者への侵襲性を抑えるため、たいていの場合、造影剤は静脈から体内へ注入する。この場合、造影剤は体内を循環した後、解析対象組織に流入する。入力関数をＣｉｎ（ｔ）、循環経路における伝達関数をｆ１ 、対象組織における伝達関数をｆ２ 、対象組織からの出力関数をＣｏｕｔ（ｔ） とおくと、Ｃｏｕｔ（ｔ） は、Ｃｉｎ（ｔ）とｆ１ を畳み込み積分したものを、更にｆ２ と畳み込み積分したものになるので、次式、
【００７３】
【数８】
【００７４】
Ｃｉｎ（ｔ）とｆ１ のたたみこみ積分は流入動脈における時間−濃度曲線（動脈入力関数）ＡＩＦ（ｔ）に、Ｃｏｕｔ（ｔ） は組織における時間−濃度曲線Ｃｒ（ｔ） に置き換えられ、次式、
【００７５】
【数９】
【００７６】
（９）式は、前述した（５）式と等価である。よって伝達関数ｆ２ は、次式、
【００７７】
【数１０】
ｆ２（ｔ） ＝Ｆ−１［Ｆ｛ＡＩＦ−１（ｔ） ｝・Ｆ｛Ｃｒ（ｔ） ｝］ …（１０）
のように表される。
【００７８】
次に、逆フィルタ演算について説明する。
【００７９】
図９は本実施の形態における逆フィルタ演算の処理手順を示すフローチャートである。同図に示すように、ＣＰＵ５により、流入動脈の時間−濃度曲線をフーリエ変換する（ステップ９０１）。次に、ＣＰＵ５により、流入動脈の時間−濃度曲線のフーリエ変換から逆フィルタを算出する（ステップ９０２）。流入動脈の時間−濃度曲線のフーリエ変換の実部をＦ｛ＡＩＦＲ｝、虚部をＦ｛ＡＩＦＩ｝とすると、逆フィルタは、次式、
【００８０】
【数１１】
【００８１】
次に、逆フィルタの高周波数成分を制御する（ステップ９０３）。ウィナーフィルタやバターワースフィルタなどの高周波数領域制御フィルタによって逆フィルタをフィルタリングすることにより、逆フィルタの高周波数成分が制御される。
【００８２】
続いて、ＣＰＵ５により、各組織の時間−濃度曲線をフーリエ変換する（ステップ９０４）。このフーリエ変換に対して、ステップ９０２で求めた逆フィルタでフィルタリングする（ステップ９０５）。このフィルタリングされた組織の時間−濃度曲線のフーリエ変換を逆フーリエ変換する伝達関数が求められる（ステップ９０６）。
【００８３】
次に、ＣＰＵ５により、ステップ９０６で求められた伝達関数から、生体機能情報を算出する（ステップ９０７）。生体機能情報の算出方法については後述する。尚、ステップ９０４からステップ９０７までは、解析する組織の断層像上での画素数分だけ繰り返す。また、ステップ９０１からステップ９０３までと、ステップ９０４は、処理の順序は逆になってもよい。
【００８４】
次に、伝達関数から生体機能情報を算出する方法について説明する。
【００８５】
図１０は生体機能情報の算出方法の説明図である。伝達関数をｆ（ｔ）とすると、ｆ（ｔ）は、例えば図１０のような形状になり、血流量（ｂｌｏｏｄｆｌｏｗ：以下、ＢＦと記述する）は、次式、
【００８６】
【数１２】
ＢＦ＝ｆｍａｘ …（１２）
のようにｆ（ｔ）の最大値から求められ、血液量（ｂｌｏｏｄｖｏｌｕｍｅ：以下、ＢＶと記述する）は、次式、
【００８７】
【数１３】
ＢＶ＝∫ｆ（ｔ）ｄｔ …（１３）
のようにｆ（ｔ）の曲線下面積から求められ、平均通過時間（ｍｅａｎｔｒａｎｓｉｔ ｔｉｍｅ： 以下、ＭＴＴ と記述する）は、次式、
【００８８】
【数１４】
ＭＴＴ ＝∫ｆ（ｔ）ｄｔ／ｆｍａｘ …（１４）
のようにｆ（ｔ）の幅から求められる。
【００８９】
ここで、ｆ ｍａｘ は、ｆ（ｔ）の最大値を表す。
【００９０】
【発明の効果】
以上説明したように本発明によれば、コンピュータ断層診断装置の提供するダイナミック撮影した断層像から生体機能情報解析を行う場合、ガンマ関数解析法や最大勾配法よりも低い造影レートで、ｄｅｃｏｎｖｏｌｕｔｉｏｎ 法よりも短い演算時間で、生体機能情報解析を行えるという効果がある。また、ＰＶＡ効果による定量性の悪化を防止できるという効果もある。更に、操作者への負担を増加させることなく不要領域の影響による無用な演算時間の増加や解析精度の悪化を防止する事ができるという効果もある。
【図面の簡単な説明】
【図１】本発明に係る血流動態解析装置のハードウェア構成例を示すブロック図
【図２】本発明に係る血流動態解析装置の断層像の入力から機能画像の表示までを示すフローチャート
【図３】不要領域を削除する方法を示すフローチャート
【図４】本実施の形態におけるＰＶＡ効果の補正方法を示す図
【図５】本実施の形態におけるＰＶＡ効果の他の補正方法を示す図
【図６】本実施の形態における血管径とピーク値の算出手順を示すフローチャート
【図７】本実施の形態における逆フィルタ法を説明するための概念図
【図８】本実施の形態における逆フィルタ法を説明するための概念図
【図９】本実施の形態における逆フィルタ演算の処理手順を示すフローチャート
【図１０】生体機能情報の算出方法の説明図
【符号の説明】
１…断層像入力装置、２…演算装置、３…インターフェイス（Ｉ／Ｆ）、４…モリ、５…中央処理装置（ＣＰＵ）、６…記録媒体、７…外部入力装置、８…表示装置、９…共通バス
Claims (3)
- コンピュータ断層像撮影手段でダイナミック撮影した断層像の入力手段と、
前記入力手段から入力した断層像に基づいて各画素ごとに画素値の時間変化の情報を表す時間−濃度曲線を求める演算手段であって、流入動脈の時間−濃度曲線及び各組織における時間−濃度曲線を求める演算手段と、
前記求めた流入動脈の時間−濃度曲線の逆フィルタを算出する手段と、
前記算出した逆フィルタを前記求めた各組織における時間−濃度曲線に掛け合わせることで各組織の伝達関数を求める手段と、
各組織の伝達関数を使用して生体機能情報を算出する手段と、
を備えたことを特徴とする血流動態解析装置。 - 前記演算手段は、断層像内の最大連結画素領域を抽出する手段と、前記抽出した最大連結画素領域を示すマスク画像を作成する手段と、前記マスク画像を使用して前記入力手段から入力する断層像のうちの当該マスク画像の領域以外の画素を削除する手段とを有し、ルームエアーや寝台、骨等の生体機能情報解析には不要な領域を自動的に除去することを特徴とする請求項１の血流動態解析装置。
- 前記演算手段は、流出静脈の時間−濃度曲線を求めるとともに、前記求めた流入動脈の時間−濃度曲線及び流出静脈の時間−濃度曲線のピーク値をそれぞれ求める手段と、前記流入動脈の時間−濃度曲線のピーク値が前記流出静脈の時間−濃度曲線のピーク値が一致するように前記流入動脈の時間−濃度曲線を補正する補正手段とを有し、流入動脈の時間−濃度曲線におけるパーシャルボリュームアベレージング効果を補正することを特徴とする請求項１又は２の血流動態解析装置。
Priority Applications (5)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2002266475A JP4013128B2 (ja) | 2002-09-12 | 2002-09-12 | 血流動態解析装置、方法、及び画像診断装置 |
| US10/527,226 US7436990B2 (en) | 2002-09-12 | 2003-09-12 | Blood flow dynamic analyzer and its method, and image diagnostic apparatus |
| PCT/JP2003/011700 WO2004024001A1 (ja) | 2002-09-12 | 2003-09-12 | 血流動態解析装置とその方法、画像診断装置 |
| CNB038215926A CN100377691C (zh) | 2002-09-12 | 2003-09-12 | 血流动态分析设备及方法和图像诊断设备 |
| EP03795424A EP1537824A4 (en) | 2002-09-12 | 2003-09-12 | DYNAMIC BLOOD FLOW ANALYZER AND ITS METHOD AND PICTOR DIAGNOSTIC |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2002266475A JP4013128B2 (ja) | 2002-09-12 | 2002-09-12 | 血流動態解析装置、方法、及び画像診断装置 |
Publications (3)
| Publication Number | Publication Date |
|---|---|
| JP2004097665A true JP2004097665A (ja) | 2004-04-02 |
| JP2004097665A5 JP2004097665A5 (ja) | 2005-11-04 |
| JP4013128B2 JP4013128B2 (ja) | 2007-11-28 |
Family
ID=31986643
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2002266475A Expired - Lifetime JP4013128B2 (ja) | 2002-09-12 | 2002-09-12 | 血流動態解析装置、方法、及び画像診断装置 |
Country Status (5)
| Country | Link |
|---|---|
| US (1) | US7436990B2 (ja) |
| EP (1) | EP1537824A4 (ja) |
| JP (1) | JP4013128B2 (ja) |
| CN (1) | CN100377691C (ja) |
| WO (1) | WO2004024001A1 (ja) |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2006025868A (ja) * | 2004-07-12 | 2006-02-02 | Ge Medical Systems Global Technology Co Llc | 画像処理装置及び画像処理方法並びにｘ線ｃｔシステム |
| JP2006263282A (ja) * | 2005-03-25 | 2006-10-05 | Hiroshima Industrial Promotion Organization | 子宮の蠕動運動を可視化するための画像処理装置、画像処理方法および画像診断支援システム、ならびにそのためのコンピュータプログラムおよびそれを記録した記録媒体 |
| JP2006326078A (ja) * | 2005-05-27 | 2006-12-07 | Hitachi Medical Corp | 血流動態解析装置、ｘ線ｃｔ装置、ｍｒｉ装置、及び血流動態解析プログラム |
| JP2007144139A (ja) * | 2005-11-02 | 2007-06-14 | Toshiba Corp | Ｘ線コンピュータ断層撮影装置及び画像処理装置 |
| JP2007536048A (ja) * | 2004-05-04 | 2007-12-13 | スティフテルセン ウニヴァズィテッフォースクニング ベルゲン | Ｍｒ画像化方法 |
| JP2007330413A (ja) * | 2006-06-14 | 2007-12-27 | Hitachi Medical Corp | 画像診断装置 |
| JP2009507545A (ja) * | 2005-09-13 | 2009-02-26 | コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ | イメージングに対する多重造影剤注入 |
| JP2009247773A (ja) * | 2008-04-10 | 2009-10-29 | Hitachi Medical Corp | 磁気共鳴イメージング装置 |
| JP2014000455A (ja) * | 2007-11-20 | 2014-01-09 | Olea Medical Sas | 患者から取得された生体画像の多数の系列を処理する方法および記録媒体 |
| JP2016101364A (ja) * | 2014-11-28 | 2016-06-02 | 株式会社東芝 | Ｘ線診断装置 |
| JP2021183165A (ja) * | 2014-03-11 | 2021-12-02 | ザ・ジョンズ・ホプキンス・ユニバーシティ | 動脈網における流量および圧力勾配を患者特定コンピュータ断層撮影アルゴリズムに基づくコントラスト分布から判断するための方法 |
Families Citing this family (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2005114229A2 (en) * | 2004-05-14 | 2005-12-01 | The General Hospital Corporation | Perfusion weighted mri with local arterial input functions |
| CN101304690B (zh) * | 2005-11-10 | 2011-06-08 | 伯拉考国际股份公司 | 基于流体力学分析在医学成像应用中检测固定造影剂 |
| DE102008020814A1 (de) | 2008-04-01 | 2009-10-15 | Siemens Aktiengesellschaft | Verfahren zur Ermittlung und Darstellung von Perfusionsparametern in der Tomographie |
| US20130039553A1 (en) * | 2010-03-09 | 2013-02-14 | Aarhus Universitet | Method for obtaining a blood flow parameter |
| GB201109344D0 (en) * | 2011-06-03 | 2011-07-20 | Siemens Medical Solutions | BIF from WB dynamic |
| EP2900142B1 (en) * | 2012-09-25 | 2019-09-04 | The Johns Hopkins University | A method for estimating flow rates, pressure gradients, coronary flow reserve, and fractional flow reserve from patient specific computed tomography angiogram-based contrast distribution data |
| CN104207797B (zh) * | 2013-06-03 | 2018-06-19 | 上海联影医疗科技有限公司 | 基于ct定位片确定病人尺寸、位置及x射线经过病人衰减的方法 |
| CN105559810B (zh) * | 2015-12-10 | 2017-08-08 | 博动医学影像科技（上海）有限公司 | 血管单位时间血流量与血流速度的计算方法 |
| KR101744424B1 (ko) * | 2015-12-24 | 2017-06-08 | (의료)길의료재단 | Mra 영상을 이용하여 관류 특성 확인용 추가 영상을 얻기 위한 방법 및 시스템 |
Family Cites Families (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US4862359A (en) * | 1984-08-31 | 1989-08-29 | Bio-Logic Systems Corporation | Topographical mapping of brain functionality from neuropsychological test results |
| IL77840A (en) * | 1986-02-10 | 1989-05-15 | Elscint Ltd | Data compression system for digital imaging |
| US4751643A (en) * | 1986-08-04 | 1988-06-14 | General Electric Company | Method and apparatus for determining connected substructures within a body |
| US5150292A (en) * | 1989-10-27 | 1992-09-22 | Arch Development Corporation | Method and system for determination of instantaneous and average blood flow rates from digital angiograms |
| JPH06133962A (ja) * | 1992-10-29 | 1994-05-17 | Hitachi Medical Corp | 画像処理装置 |
| JP2001511054A (ja) * | 1997-12-08 | 2001-08-07 | コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ | 磁気共鳴により対象を画像化する方法及び装置 |
| EP1171028B1 (en) * | 1999-03-26 | 2005-11-16 | Leif Ostergaard | System for determining haemodynamic indices by use of tomographic data |
| EP2281503B1 (en) * | 1999-09-24 | 2015-04-29 | National Research Council of Canada | Method for performing intra-operative angiography |
| US6618609B2 (en) * | 2001-03-30 | 2003-09-09 | Koninklifke Philips Electronics, N.V. | Flow adaptive guided determination of imaging parameters for 3D magnetic resonance angiography |
| US6546275B2 (en) * | 2001-06-25 | 2003-04-08 | Wisconsin Alumni Research Foundation | Determination of the arterial input function in dynamic contrast-enhanced MRI |
| WO2004089218A1 (ja) * | 2003-04-04 | 2004-10-21 | Hitachi Medical Corporation | 機能画像の表示方法及び装置 |
| CA2653068A1 (en) * | 2006-07-13 | 2008-01-17 | The Regents Of The University Of Colorado | Echo particle image velocity (epiv) and echo particle tracking velocimetry (eptv) system and method |
-
2002
- 2002-09-12 JP JP2002266475A patent/JP4013128B2/ja not_active Expired - Lifetime
-
2003
- 2003-09-12 US US10/527,226 patent/US7436990B2/en active Active
- 2003-09-12 EP EP03795424A patent/EP1537824A4/en not_active Withdrawn
- 2003-09-12 CN CNB038215926A patent/CN100377691C/zh not_active Expired - Lifetime
- 2003-09-12 WO PCT/JP2003/011700 patent/WO2004024001A1/ja active Application Filing
Cited By (14)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2007536048A (ja) * | 2004-05-04 | 2007-12-13 | スティフテルセン ウニヴァズィテッフォースクニング ベルゲン | Ｍｒ画像化方法 |
| JP2006025868A (ja) * | 2004-07-12 | 2006-02-02 | Ge Medical Systems Global Technology Co Llc | 画像処理装置及び画像処理方法並びにｘ線ｃｔシステム |
| JP4535795B2 (ja) * | 2004-07-12 | 2010-09-01 | ジーイー・メディカル・システムズ・グローバル・テクノロジー・カンパニー・エルエルシー | 画像処理装置及びｘ線ｃｔシステム |
| JP2006263282A (ja) * | 2005-03-25 | 2006-10-05 | Hiroshima Industrial Promotion Organization | 子宮の蠕動運動を可視化するための画像処理装置、画像処理方法および画像診断支援システム、ならびにそのためのコンピュータプログラムおよびそれを記録した記録媒体 |
| JP4575821B2 (ja) * | 2005-03-25 | 2010-11-04 | 広島市 | 子宮の蠕動運動を可視化するための画像処理装置、画像処理方法および画像診断支援システム、ならびにそのためのコンピュータプログラムおよびそれを記録した記録媒体 |
| JP2006326078A (ja) * | 2005-05-27 | 2006-12-07 | Hitachi Medical Corp | 血流動態解析装置、ｘ線ｃｔ装置、ｍｒｉ装置、及び血流動態解析プログラム |
| JP2009507545A (ja) * | 2005-09-13 | 2009-02-26 | コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ | イメージングに対する多重造影剤注入 |
| JP2007144139A (ja) * | 2005-11-02 | 2007-06-14 | Toshiba Corp | Ｘ線コンピュータ断層撮影装置及び画像処理装置 |
| JP2007330413A (ja) * | 2006-06-14 | 2007-12-27 | Hitachi Medical Corp | 画像診断装置 |
| JP2014000455A (ja) * | 2007-11-20 | 2014-01-09 | Olea Medical Sas | 患者から取得された生体画像の多数の系列を処理する方法および記録媒体 |
| JP2009247773A (ja) * | 2008-04-10 | 2009-10-29 | Hitachi Medical Corp | 磁気共鳴イメージング装置 |
| JP2021183165A (ja) * | 2014-03-11 | 2021-12-02 | ザ・ジョンズ・ホプキンス・ユニバーシティ | 動脈網における流量および圧力勾配を患者特定コンピュータ断層撮影アルゴリズムに基づくコントラスト分布から判断するための方法 |
| JP7303260B2 (ja) | 2014-03-11 | 2023-07-04 | ザ・ジョンズ・ホプキンス・ユニバーシティ | 動脈網における流量および圧力勾配を患者特定コンピュータ断層撮影アルゴリズムに基づくコントラスト分布から判断するための方法 |
| JP2016101364A (ja) * | 2014-11-28 | 2016-06-02 | 株式会社東芝 | Ｘ線診断装置 |
Also Published As
| Publication number | Publication date |
|---|---|
| WO2004024001A1 (ja) | 2004-03-25 |
| US20060098852A1 (en) | 2006-05-11 |
| JP4013128B2 (ja) | 2007-11-28 |
| EP1537824A1 (en) | 2005-06-08 |
| US7436990B2 (en) | 2008-10-14 |
| EP1537824A4 (en) | 2009-04-15 |
| CN100377691C (zh) | 2008-04-02 |
| CN1681438A (zh) | 2005-10-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP4013128B2 (ja) | 血流動態解析装置、方法、及び画像診断装置 | |
| JP4193168B2 (ja) | 血流動態解析装置及び方法 | |
| US7676257B2 (en) | Method and apparatus for segmenting structure in CT angiography | |
| JP6559678B2 (ja) | スペクトル画像データ処理 | |
| EP2399238B1 (en) | Functional imaging | |
| RU2717885C1 (ru) | Оценка потока, сопротивления или давления на основании измерений давления или потока и ангиографии | |
| JP3656695B2 (ja) | 骨計測方法および装置 | |
| WO2010055405A2 (en) | Method and system for mapping tissue status of acute stroke | |
| JP4804039B2 (ja) | 血流動態解析装置、ｘ線ｃｔ装置、ｍｒｉ装置、及び血流動態解析プログラム | |
| US20120170820A1 (en) | Methods and apparatus for comparing 3d and 2d image data | |
| CN109716388B (zh) | 图像数据中的降噪 | |
| Li et al. | A robust noise reduction technique for time resolved CT | |
| Wu et al. | Cortical bone vessel identification and quantification on contrast-enhanced MR images | |
| Bruder et al. | Spatio-temporal filtration of dynamic CT data using diffusion filters | |
| US8290224B2 (en) | Method and device for imaging cyclically moving objects | |
| Eck et al. | Comparison of quantitative myocardial perfusion imaging CT to fluorescent microsphere-based flow from high-resolution cryo-images | |
| JP2003299646A (ja) | 画像解析装置 | |
| CN112168191A (zh) | 从第一三维计算机断层成像数据组提供分析数据组的方法 | |
| JP7191038B2 (ja) | スペクトルイメージングｆｆｒ | |
| JP5337416B2 (ja) | 画像処理装置及び画像診断装置 | |
| JP4964450B2 (ja) | 画像処理方法、画像処理プログラム、及び画像処理装置 | |
| Sun et al. | Left atrium wall-mapping application for wall thickness visualisation | |
| Elfarnawany et al. | Improved objective selection of power Doppler wall-filter cut-off velocity for accurate vascular quantification | |
| CN108289650A (zh) | 计算机断层摄影成像方法 | |
| JPH11155849A (ja) | 骨画像情報取得方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20050902 |
|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20050902 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20070816 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20070829 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100921 Year of fee payment: 3 |
|
| R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 Ref document number: 4013128 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100921 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110921 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110921 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120921 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120921 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130921 Year of fee payment: 6 |
|
| S111 | Request for change of ownership or part of ownership |
Free format text: JAPANESE INTERMEDIATE CODE: R313111 |
|
| S533 | Written request for registration of change of name |
Free format text: JAPANESE INTERMEDIATE CODE: R313533 |
|
| R350 | Written notification of registration of transfer |
Free format text: JAPANESE INTERMEDIATE CODE: R350 |
|
| S111 | Request for change of ownership or part of ownership |
Free format text: JAPANESE INTERMEDIATE CODE: R313111 |
|
| R350 | Written notification of registration of transfer |
Free format text: JAPANESE INTERMEDIATE CODE: R350 |
|
| EXPY | Cancellation because of completion of term |