JP4013128B2 - 血流動態解析装置、方法、及び画像診断装置 - Google Patents
血流動態解析装置、方法、及び画像診断装置 Download PDFInfo
- Publication number
- JP4013128B2 JP4013128B2 JP2002266475A JP2002266475A JP4013128B2 JP 4013128 B2 JP4013128 B2 JP 4013128B2 JP 2002266475 A JP2002266475 A JP 2002266475A JP 2002266475 A JP2002266475 A JP 2002266475A JP 4013128 B2 JP4013128 B2 JP 4013128B2
- Authority
- JP
- Japan
- Prior art keywords
- time
- pixel
- tomographic image
- concentration curve
- blood flow
- 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.)
- Expired - Lifetime
Links
- 238000004458 analytical method Methods 0.000 title claims description 51
- 238000000034 method Methods 0.000 title claims description 46
- 230000017531 blood circulation Effects 0.000 title claims description 44
- 238000002059 diagnostic imaging Methods 0.000 title claims description 4
- 210000001367 artery Anatomy 0.000 claims description 64
- 238000004364 calculation method Methods 0.000 claims description 53
- 230000000694 effects Effects 0.000 claims description 32
- 210000003462 vein Anatomy 0.000 claims description 29
- 238000012546 transfer Methods 0.000 claims description 23
- 238000003384 imaging method Methods 0.000 claims description 13
- 210000000988 bone and bone Anatomy 0.000 claims description 6
- 238000012937 correction Methods 0.000 claims description 6
- 230000002123 temporal effect Effects 0.000 claims description 6
- 238000012935 Averaging Methods 0.000 claims description 5
- 238000000605 extraction Methods 0.000 claims description 5
- 230000000004 hemodynamic effect Effects 0.000 claims description 3
- 238000003325 tomography Methods 0.000 claims 1
- 230000006870 function Effects 0.000 description 47
- 210000004204 blood vessel Anatomy 0.000 description 34
- 230000008827 biological function Effects 0.000 description 21
- 210000001519 tissue Anatomy 0.000 description 20
- 238000010586 diagram Methods 0.000 description 12
- 238000012545 processing Methods 0.000 description 11
- 238000002591 computed tomography Methods 0.000 description 9
- 238000012905 input function Methods 0.000 description 9
- 239000002872 contrast media Substances 0.000 description 8
- 238000003745 diagnosis Methods 0.000 description 7
- 210000000056 organ Anatomy 0.000 description 7
- 230000006866 deterioration Effects 0.000 description 4
- 239000000284 extract Substances 0.000 description 4
- 239000008280 blood Substances 0.000 description 3
- 210000004369 blood Anatomy 0.000 description 3
- 230000004087 circulation Effects 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 238000011002 quantification Methods 0.000 description 3
- 230000001154 acute effect Effects 0.000 description 2
- 210000004556 brain 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
- 238000010191 image analysis Methods 0.000 description 2
- 238000002372 labelling Methods 0.000 description 2
- 238000009206 nuclear medicine Methods 0.000 description 2
- 230000010412 perfusion Effects 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
- 210000002551 anterior cerebral artery Anatomy 0.000 description 1
- 230000002490 cerebral effect Effects 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 238000002405 diagnostic procedure Methods 0.000 description 1
- 238000013213 extrapolation Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- UIWYJDYFSGRHKR-UHFFFAOYSA-N gadolinium atom Chemical compound [Gd] UIWYJDYFSGRHKR-UHFFFAOYSA-N 0.000 description 1
- 238000001802 infusion Methods 0.000 description 1
- 229910052740 iodine Inorganic materials 0.000 description 1
- 239000011630 iodine Substances 0.000 description 1
- 230000003902 lesion Effects 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 210000003657 middle cerebral artery Anatomy 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
- 238000000926 separation method Methods 0.000 description 1
- 210000000798 superior sagittal sinus Anatomy 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)
- Surgery (AREA)
- Public Health (AREA)
- Biophysics (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Veterinary Medicine (AREA)
- Pathology (AREA)
- Radiology & Medical Imaging (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Molecular Biology (AREA)
- Animal Behavior & Ethology (AREA)
- High Energy & Nuclear Physics (AREA)
- Optics & Photonics (AREA)
- Oral & Maxillofacial Surgery (AREA)
- Dentistry (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Quality & Reliability (AREA)
- Vascular Medicine (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Hematology (AREA)
- Cardiology (AREA)
- Physiology (AREA)
- Apparatus For Radiation Diagnosis (AREA)
- Image Processing (AREA)
- Magnetic Resonance Imaging Apparatus (AREA)
Description
【０００１】
【発明の属する技術分野】
本発明は血流動態解析装置に係り、特にＸ線や電磁波を用いて被検体の時間変化に依存する生体機能情報を得ることが可能なコンピュータ断層診断装置の提供する断層像から生体機能情報解析等の画像解析を行うことに対して有用な技術に関する。
【０００２】
【従来の技術】
従来、ダイナミック撮影を行うコンピュータ断層診断装置として、陽電子放射断層撮影装置（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)とｆを畳み込み積分したものになり、次式、
【００６１】
【数４】
のように表される。
【００６２】
造影剤を動脈から体内に注入する場合、Cin(t)は流入動脈における時間−濃度曲線（動脈入力関数）AIF(t)に、Cout(t) は組織における時間−濃度曲線Cr(t) に置き換えられ、次式、
【００６３】
【数５】
のように表される。
【００６４】
（５）式を周波数空間に変換すると、次式、
【００６５】
【数６】
Ｆ｛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 と畳み込み積分したものになるので、次式、
【００７３】
【数８】
のようになる。
【００７４】
Cin(t)とｆ1 のたたみこみ積分は流入動脈における時間−濃度曲線（動脈入力関数）AIF(t)に、Cout(t) は組織における時間−濃度曲線Cr(t) に置き換えられ、次式、
【００７５】
【数９】
のように表される。
【００７６】
（９）式は、前述した（５）式と等価である。よって伝達関数ｆ2 は、次式、
【００７７】
【数１０】
f2(t) ＝Ｆ-1［Ｆ｛AIF-1(t) ｝・Ｆ｛Cr(t) ｝］ …（１０）
のように表される。
【００７８】
次に、逆フィルタ演算について説明する。
【００７９】
図９は本実施の形態における逆フィルタ演算の処理手順を示すフローチャートである。同図に示すように、ＣＰＵ５により、流入動脈の時間−濃度曲線をフーリエ変換する（ステップ９０１）。次に、ＣＰＵ５により、流入動脈の時間−濃度曲線のフーリエ変換から逆フィルタを算出する（ステップ９０２）。流入動脈の時間−濃度曲線のフーリエ変換の実部をＦ｛AIFR｝、虚部をＦ｛AIFI｝とすると、逆フィルタは、次式、
【００８０】
【数１１】
で表される。
【００８１】
次に、逆フィルタの高周波数成分を制御する（ステップ９０３）。ウィナーフィルタやバターワースフィルタなどの高周波数領域制御フィルタによって逆フィルタをフィルタリングすることにより、逆フィルタの高周波数成分が制御される。
【００８２】
続いて、ＣＰＵ５により、各組織の時間−濃度曲線をフーリエ変換する（ステップ９０４）。このフーリエ変換に対して、ステップ９０２で求めた逆フィルタでフィルタリングする（ステップ９０５）。このフィルタリングされた組織の時間−濃度曲線のフーリエ変換を逆フーリエ変換する伝達関数が求められる（ステップ９０６）。
【００８３】
次に、ＣＰＵ５により、ステップ９０６で求められた伝達関数から、生体機能情報を算出する（ステップ９０７）。生体機能情報の算出方法については後述する。尚、ステップ９０４からステップ９０７までは、解析する組織の断層像上での画素数分だけ繰り返す。また、ステップ９０１からステップ９０３までと、ステップ９０４は、処理の順序は逆になってもよい。
【００８４】
次に、伝達関数から生体機能情報を算出する方法について説明する。
【００８５】
図１０は生体機能情報の算出方法の説明図である。伝達関数を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 法よりも短い演算時間で、生体機能情報解析を行えるという効果がある。また、ＰＶＡ効果による定量性の悪化を防止できるという効果もある。更に、操作者への負担を増加させることなく不要領域の影響による無用な演算時間の増加や解析精度の悪化を防止する事ができるという効果もある。
【図面の簡単な説明】
【図１】 本発明に係る血流動態解析装置のハードウェア構成例を示すブロック図
【図２】 本発明に係る血流動態解析装置の断層像の入力から機能画像の表示までを示すフローチャート
【図３】 不要領域を削除する方法を示すフローチャート
【図４】 本実施の形態におけるＰＶＡ効果の補正方法を示す図
【図５】 本実施の形態におけるＰＶＡ効果の他の補正方法を示す図
【図６】 本実施の形態における血管径とピーク値の算出手順を示すフローチャート
【図７】 本実施の形態における逆フィルタ法を説明するための概念図
【図８】 本実施の形態における逆フィルタ法を説明するための概念図
【図９】 本実施の形態における逆フィルタ演算の処理手順を示すフローチャート
【図１０】 生体機能情報の算出方法の説明図
【符号の説明】
１…断層像入力装置、２…演算装置、３…インターフェイス（Ｉ／Ｆ）、４…モリ、５…中央処理装置（ＣＰＵ）、６…記録媒体、７…外部入力装置、８…表示装置、９…共通バス
【発明の属する技術分野】
本発明は血流動態解析装置に係り、特にＸ線や電磁波を用いて被検体の時間変化に依存する生体機能情報を得ることが可能なコンピュータ断層診断装置の提供する断層像から生体機能情報解析等の画像解析を行うことに対して有用な技術に関する。
【０００２】
【従来の技術】
従来、ダイナミック撮影を行うコンピュータ断層診断装置として、陽電子放射断層撮影装置（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)とｆを畳み込み積分したものになり、次式、
【００６１】
【数４】
【００６２】
造影剤を動脈から体内に注入する場合、Cin(t)は流入動脈における時間−濃度曲線（動脈入力関数）AIF(t)に、Cout(t) は組織における時間−濃度曲線Cr(t) に置き換えられ、次式、
【００６３】
【数５】
【００６４】
（５）式を周波数空間に変換すると、次式、
【００６５】
【数６】
Ｆ｛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 と畳み込み積分したものになるので、次式、
【００７３】
【数８】
【００７４】
Cin(t)とｆ1 のたたみこみ積分は流入動脈における時間−濃度曲線（動脈入力関数）AIF(t)に、Cout(t) は組織における時間−濃度曲線Cr(t) に置き換えられ、次式、
【００７５】
【数９】
【００７６】
（９）式は、前述した（５）式と等価である。よって伝達関数ｆ2 は、次式、
【００７７】
【数１０】
f2(t) ＝Ｆ-1［Ｆ｛AIF-1(t) ｝・Ｆ｛Cr(t) ｝］ …（１０）
のように表される。
【００７８】
次に、逆フィルタ演算について説明する。
【００７９】
図９は本実施の形態における逆フィルタ演算の処理手順を示すフローチャートである。同図に示すように、ＣＰＵ５により、流入動脈の時間−濃度曲線をフーリエ変換する（ステップ９０１）。次に、ＣＰＵ５により、流入動脈の時間−濃度曲線のフーリエ変換から逆フィルタを算出する（ステップ９０２）。流入動脈の時間−濃度曲線のフーリエ変換の実部をＦ｛AIFR｝、虚部をＦ｛AIFI｝とすると、逆フィルタは、次式、
【００８０】
【数１１】
【００８１】
次に、逆フィルタの高周波数成分を制御する（ステップ９０３）。ウィナーフィルタやバターワースフィルタなどの高周波数領域制御フィルタによって逆フィルタをフィルタリングすることにより、逆フィルタの高周波数成分が制御される。
【００８２】
続いて、ＣＰＵ５により、各組織の時間−濃度曲線をフーリエ変換する（ステップ９０４）。このフーリエ変換に対して、ステップ９０２で求めた逆フィルタでフィルタリングする（ステップ９０５）。このフィルタリングされた組織の時間−濃度曲線のフーリエ変換を逆フーリエ変換する伝達関数が求められる（ステップ９０６）。
【００８３】
次に、ＣＰＵ５により、ステップ９０６で求められた伝達関数から、生体機能情報を算出する（ステップ９０７）。生体機能情報の算出方法については後述する。尚、ステップ９０４からステップ９０７までは、解析する組織の断層像上での画素数分だけ繰り返す。また、ステップ９０１からステップ９０３までと、ステップ９０４は、処理の順序は逆になってもよい。
【００８４】
次に、伝達関数から生体機能情報を算出する方法について説明する。
【００８５】
図１０は生体機能情報の算出方法の説明図である。伝達関数を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 法よりも短い演算時間で、生体機能情報解析を行えるという効果がある。また、ＰＶＡ効果による定量性の悪化を防止できるという効果もある。更に、操作者への負担を増加させることなく不要領域の影響による無用な演算時間の増加や解析精度の悪化を防止する事ができるという効果もある。
【図面の簡単な説明】
【図１】 本発明に係る血流動態解析装置のハードウェア構成例を示すブロック図
【図２】 本発明に係る血流動態解析装置の断層像の入力から機能画像の表示までを示すフローチャート
【図３】 不要領域を削除する方法を示すフローチャート
【図４】 本実施の形態におけるＰＶＡ効果の補正方法を示す図
【図５】 本実施の形態におけるＰＶＡ効果の他の補正方法を示す図
【図６】 本実施の形態における血管径とピーク値の算出手順を示すフローチャート
【図７】 本実施の形態における逆フィルタ法を説明するための概念図
【図８】 本実施の形態における逆フィルタ法を説明するための概念図
【図９】 本実施の形態における逆フィルタ演算の処理手順を示すフローチャート
【図１０】 生体機能情報の算出方法の説明図
【符号の説明】
１…断層像入力装置、２…演算装置、３…インターフェイス（Ｉ／Ｆ）、４…モリ、５…中央処理装置（ＣＰＵ）、６…記録媒体、７…外部入力装置、８…表示装置、９…共通バス
Claims (15)
- 断層像撮影手段により撮影された断層像を入力する入力手段と、
この入力手段により入力された断層像の各画素の時間変化の情報を表す時間−濃度曲線を求める第１の演算手段と、
この第１の演算手段により求められた断層像の各画素の時間−濃度曲線から流入動脈の時間−濃度曲線を抽出する第２の演算手段と、
この第２の演算手段により抽出された流入動脈の時間−濃度曲線から逆フィルタ関数を求める第３の演算手段と、
この第３の演算手段により求められた逆フィルタ関数と前記第１の演算手段により求められた断層像の各画素の時間−濃度曲線に基づき前記断層像の各画素の伝達関数を求める第４の演算手段と、
この第４の演算手段により求められた前記断層像の各画素の伝達関数を用いて血流動態解析像を求める第５の演算手段と、
を備えたことを特徴とする血流動態解析装置。 - 前記第３の演算手段は、前記求められた逆フィルタ関数の高周波成分を制御する高周波制御フィルタ手段を備えたことを特徴とする請求項１に記載の血流動態解析装置。
- 前記高周波制御フィルタ手段は、前記第２の演算手段に抽出された流入動脈の時間―濃度曲線からフィルタ関数を求めるフィルタ関数演算手段を備えたことを特徴とする請求項２に記載の血流動態解析装置。
- 前記入力手段により入力された断層像のうちの最大連結画素領域を抽出する手段と、この抽出手段により抽出された最大連結画素領域に基づき前記入力手段から入力される断層像のうち血流動態解析に不要な領域を除去する手段と、を備えたことを特徴とする請求項１に記載の血流動態解析装置。
- 前記不要領域は、ルームエアー、寝台、及び骨を含むことを特徴とする請求項４に記載の血流動態解析装置。
- 前記入力手段により入力された断層像の各画素の時間−濃度曲線から流出静脈の時間−濃度曲線を求めるとともに、その求められた流入動脈の時間−濃度曲線及び前記流出静脈の時間−濃度曲線のピーク値をそれぞれ求める手段と、前記求められた流入動脈の時間−濃度曲線のピーク値と前記流出静脈の時間―濃度曲線のピーク値とが一致するように流入動脈の時間−濃度曲線におけるパーシャルボリュームアベレージング効果を補正する補正手段と、を備えたことを特徴とする請求項１に記載の血流動態解析装置。
- 断層像撮影装置により撮影された断層像を入力する入力ステップと、
この入力ステップにより入力された断層像の各画素の時間変化の情報を表す時間−濃度曲線を求める第１の演算ステップと、
この第１の演算ステップにより求められた断層像の各画素の時間−濃度曲線から流入動脈の時間−濃度曲線を抽出する第２の演算ステップと、
この第２の演算ステップにより抽出された流入動脈の時間−濃度曲線から逆フィルタ関数を求める第３の演算ステップと、
この第３の演算ステップにより求められた逆フィルタ関数と前記第１の演算ステップにより求められた断層像の各画素の時間−濃度曲線に基づき前記断層像の各画素の伝達関数を求める第４の演算ステップと、
この第４の演算ステップにより求められた前記断層像の各画素の伝達関数を用いて血流動態解析像を求める第５の演算ステップと、
を含むことを特徴とする血流動態解析方法。 - 前記第３の演算ステップは、前記求められた逆フィルタ関数の高周波成分を制御する高周波制御フィルタステップを含むことを特徴とする請求項７に記載の血流動態解析方法。
- 前記高周波制御フィルタステップは、前記第２の演算ステップに抽出された流入動脈の時間―濃度曲線からフィルタ関数を求めるフィルタ関数演算ステップを含むことを特徴とする請求項８に記載の血流動態解析方法。
- 前記入力ステップにより入力された断層像のうちの最大連結画素領域を抽出するステップと、この抽出ステップにより抽出された最大連結画素領域に基づき前記入力ステップから入力される断層像のうち血流動態解析に不要な領域を除去するステップと、を含むことを特徴とする請求項７に記載の血流動態解析方法。
- 前記不要領域は、ルームエアー、寝台、及び骨を含むことを特徴とする請求項１０に記載の血流動態解析方法。
- 前記入力ステップにより入力された断層像の各画素の時間−濃度曲線から流出静脈の時間−濃度曲線を求めるとともに、その求められた流入動脈の時間−濃度曲線及び前記流出静脈の時間−濃度曲線のピーク値をそれぞれ求めるステップと、前記求められた流入動脈の時間−濃度曲線のピーク値と前記流出静脈の時間―濃度曲線のピーク値とが一致するように流入動脈の時間−濃度曲線におけるパーシャルボリュームアベレージング効果を補正する補正ステップと、を含むことを特徴とする請求項７に記載の血流動態解析方法。
- 生体から断層像を撮影可能な画像診断装置において、
前記撮影された断層像を入力する入力手段と、
この入力手段により入力された断層像の各画素の時間変化の情報を表す時間−濃度曲線を求める第１の演算手段と、
この第１の演算手段により求められた断層像の各画素の時間−濃度曲線から流入動脈の時間−濃度曲線を抽出する第２の演算手段と、
この第２の演算手段により抽出された流入動脈の時間−濃度曲線から逆フィルタ関数を求める第３の演算手段と、
この第３の演算手段により求められた逆フィルタ関数と前記第１の演算手段により求められた断層像の各画素の時間−濃度曲線に基づき前記断層像の各画素の伝達関数を求める第４の演算手段と、
この第４の演算手段により求められた前記断層像の各画素の伝達関数を用いて血流動態解析像を求める第５の演算手段と、
を備えたことを特徴とする画像診断装置。 - 前記第３の演算手段は、前記求められた逆フィルタ関数の高周波成分を制御する高周波制御フィルタ手段を備えたことを特徴とする請求項１３に記載の画像診断装置。
- 前記高周波制御フィルタ手段は、前記第２の演算手段に抽出された流入動脈の時間―濃度曲線からフィルタ関数を求めるフィルタ関数演算手段を備えたことを特徴とする請求項１４に記載の画像診断装置。
Priority Applications (5)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2002266475A JP4013128B2 (ja) | 2002-09-12 | 2002-09-12 | 血流動態解析装置、方法、及び画像診断装置 |
| PCT/JP2003/011700 WO2004024001A1 (ja) | 2002-09-12 | 2003-09-12 | 血流動態解析装置とその方法、画像診断装置 |
| US10/527,226 US7436990B2 (en) | 2002-09-12 | 2003-09-12 | Blood flow dynamic analyzer and its method, and image diagnostic apparatus |
| EP03795424A EP1537824A4 (en) | 2002-09-12 | 2003-09-12 | DYNAMIC BLOOD FLOW ANALYZER AND ITS METHOD AND PICTOR DIAGNOSTIC |
| CNB038215926A CN100377691C (zh) | 2002-09-12 | 2003-09-12 | 血流动态分析设备及方法和图像诊断设备 |
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 JP2004097665A (ja) | 2004-04-02 |
| JP2004097665A5 JP2004097665A5 (ja) | 2005-11-04 |
| JP4013128B2 true 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) |
Families Citing this family (20)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101002104B (zh) * | 2004-05-04 | 2012-01-18 | 卑尔根大学研究基金会 | 动脉输入和组织残留函数在灌注磁共振成像中的盲确定 |
| US7885442B2 (en) * | 2004-05-14 | 2011-02-08 | The General Hospital Corporation | Perfusion weighted MRI with local arterial input functions |
| JP4535795B2 (ja) * | 2004-07-12 | 2010-09-01 | ジーイー・メディカル・システムズ・グローバル・テクノロジー・カンパニー・エルエルシー | 画像処理装置及びｘ線ｃｔシステム |
| JP4575821B2 (ja) * | 2005-03-25 | 2010-11-04 | 広島市 | 子宮の蠕動運動を可視化するための画像処理装置、画像処理方法および画像診断支援システム、ならびにそのためのコンピュータプログラムおよびそれを記録した記録媒体 |
| JP4804039B2 (ja) * | 2005-05-27 | 2011-10-26 | 株式会社日立メディコ | 血流動態解析装置、ｘ線ｃｔ装置、ｍｒｉ装置、及び血流動態解析プログラム |
| EP1926425B1 (en) | 2005-09-13 | 2015-08-19 | Philips Intellectual Property & Standards GmbH | Multiple contrast agent injection for imaging |
| JP2007144139A (ja) * | 2005-11-02 | 2007-06-14 | Toshiba Corp | Ｘ線コンピュータ断層撮影装置及び画像処理装置 |
| WO2007054544A1 (en) * | 2005-11-10 | 2007-05-18 | Bracco Research Sa | Detection of immobilized contrast agent in medical imaging applications based on flow dynamics analysis |
| JP2007330413A (ja) * | 2006-06-14 | 2007-12-27 | Hitachi Medical Corp | 画像診断装置 |
| US8320647B2 (en) * | 2007-11-20 | 2012-11-27 | Olea Medical | Method and system for processing multiple series of biological images obtained from a patient |
| DE102008020814A1 (de) | 2008-04-01 | 2009-10-15 | Siemens Aktiengesellschaft | Verfahren zur Ermittlung und Darstellung von Perfusionsparametern in der Tomographie |
| JP5352109B2 (ja) * | 2008-04-10 | 2013-11-27 | 株式会社日立メディコ | 磁気共鳴イメージング装置 |
| WO2011110174A1 (en) * | 2010-03-09 | 2011-09-15 | Aarhus Universitet | A 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射线经过病人衰减的方法 |
| WO2015138555A2 (en) * | 2014-03-11 | 2015-09-17 | The Johns Hopkins University | A method for estimating flow rates and pressure gradients in arterial networks from patient specific computed tomography angiogram-based contrast distribution data |
| JP6419551B2 (ja) * | 2014-11-28 | 2018-11-07 | キヤノンメディカルシステムズ株式会社 | Ｘ線診断装置、画像処理装置及び画像処理プログラム |
| 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 | 画像処理装置 |
| DE69841476D1 (de) * | 1997-12-08 | 2010-03-18 | Koninkl Philips Electronics Nv | Gerät und verfahren zur bildformung eines objektes durch magnetische resonanz |
| WO2000057777A1 (en) * | 1999-03-26 | 2000-10-05 | Oestergaard Leif | Method for determining haemodynamic indices by use of tomographic data |
| JP3881550B2 (ja) * | 1999-09-24 | 2007-02-14 | ナショナル・リサーチ・カウンシル・オブ・カナダ | 術中血管造影を行なうための装置 |
| 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 |
| JP4558645B2 (ja) * | 2003-04-04 | 2010-10-06 | 株式会社日立メディコ | 画像表示方法及び装置 |
| JP2010503421A (ja) * | 2006-07-13 | 2010-02-04 | ザ リージェンツ オブ ザ ユニバーシティ オブ コロラド | エコー粒子画像速度（ｅｐｉｖ）およびエコー粒子追跡速度測定（ｅｐｔｖ）システムおよび方法 |
-
2002
- 2002-09-12 JP JP2002266475A patent/JP4013128B2/ja not_active Expired - Lifetime
-
2003
- 2003-09-12 EP EP03795424A patent/EP1537824A4/en not_active Withdrawn
- 2003-09-12 US US10/527,226 patent/US7436990B2/en active Active
- 2003-09-12 WO PCT/JP2003/011700 patent/WO2004024001A1/ja active Application Filing
- 2003-09-12 CN CNB038215926A patent/CN100377691C/zh not_active Expired - Lifetime
Also Published As
| Publication number | Publication date |
|---|---|
| EP1537824A1 (en) | 2005-06-08 |
| JP2004097665A (ja) | 2004-04-02 |
| CN1681438A (zh) | 2005-10-12 |
| US20060098852A1 (en) | 2006-05-11 |
| US7436990B2 (en) | 2008-10-14 |
| EP1537824A4 (en) | 2009-04-15 |
| WO2004024001A1 (ja) | 2004-03-25 |
| CN100377691C (zh) | 2008-04-02 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP4013128B2 (ja) | 血流動態解析装置、方法、及び画像診断装置 | |
| JP4193168B2 (ja) | 血流動態解析装置及び方法 | |
| JP6559678B2 (ja) | スペクトル画像データ処理 | |
| EP2399238B1 (en) | Functional imaging | |
| EP2375969B1 (en) | Method and system for mapping tissue status of acute stroke | |
| RU2717885C1 (ru) | Оценка потока, сопротивления или давления на основании измерений давления или потока и ангиографии | |
| US20050113679A1 (en) | Method and apparatus for segmenting structure in CT angiography | |
| Erlandsson et al. | MR imaging–guided partial volume correction of PET data in PET/MR imaging | |
| WO2005046478A1 (ja) | 画像処理方法、画像処理装置、医用画像診断支援システム、及び時間軸方向フィルタリング方法 | |
| JP4804039B2 (ja) | 血流動態解析装置、ｘ線ｃｔ装置、ｍｒｉ装置、及び血流動態解析プログラム | |
| JP7303260B2 (ja) | 動脈網における流量および圧力勾配を患者特定コンピュータ断層撮影アルゴリズムに基づくコントラスト分布から判断するための方法 | |
| JP5801316B2 (ja) | 副側血流評価 | |
| EP3455824B1 (en) | Noise reduction in image data | |
| Herzog et al. | CT perfusion imaging of the lung in pulmonary embolism1 | |
| Fang et al. | Improving low-dose blood–brain barrier permeability quantification using sparse high-dose induced prior for patlak model | |
| Bruder et al. | Spatio-temporal filtration of dynamic CT data using diffusion filters | |
| Li et al. | A robust noise reduction technique for time resolved CT | |
| Sun et al. | Left atrium wall-mapping application for wall thickness visualisation | |
| Sun et al. | Improving the diagnostic performance of computed tomography angiography for intracranial large arterial stenosis by a novel super-resolution algorithm based on multi-scale residual denoising generative adversarial network | |
| Levi et al. | Comparison of automated beam hardening correction (ABHC) algorithms for myocardial perfusion imaging using computed tomography | |
| CN108289650A (zh) | 计算机断层摄影成像方法 | |
| Bindschadler et al. | Quantitative myocardial perfusion from static cardiac and dynamic arterial CT | |
| Murayama et al. | A Bayesian estimation method for cerebral blood flow measurement by area-detector CT perfusion imaging | |
| Sakaguchi et al. | Development of a method for automated and stable myocardial perfusion measurement using coronary X-ray angiography images | |
| Caldeira et al. | MR-based Corrections for Quantitative PET Image |
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 |