JP4193168B2 - 血流動態解析装置及び方法 - Google Patents
血流動態解析装置及び方法 Download PDFInfo
- Publication number
- JP4193168B2 JP4193168B2 JP2002025309A JP2002025309A JP4193168B2 JP 4193168 B2 JP4193168 B2 JP 4193168B2 JP 2002025309 A JP2002025309 A JP 2002025309A JP 2002025309 A JP2002025309 A JP 2002025309A JP 4193168 B2 JP4193168 B2 JP 4193168B2
- Authority
- JP
- Japan
- Prior art keywords
- pixel
- change information
- blood flow
- time change
- filter
- 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
- 230000017531 blood circulation Effects 0.000 title claims description 62
- 238000000034 method Methods 0.000 title claims description 60
- 230000008859 change Effects 0.000 claims description 45
- 239000011159 matrix material Substances 0.000 claims description 42
- 238000004458 analytical method Methods 0.000 claims description 32
- 238000009499 grossing Methods 0.000 claims description 30
- 230000008569 process Effects 0.000 claims description 21
- 238000003384 imaging method Methods 0.000 claims description 20
- 238000012545 processing Methods 0.000 claims description 19
- 238000001914 filtration Methods 0.000 claims description 7
- 230000002123 temporal effect Effects 0.000 claims description 6
- 238000002591 computed tomography Methods 0.000 claims description 5
- 238000005206 flow analysis Methods 0.000 claims description 3
- 230000006870 function Effects 0.000 description 17
- 238000004364 calculation method Methods 0.000 description 11
- 238000002595 magnetic resonance imaging Methods 0.000 description 10
- 238000004422 calculation algorithm Methods 0.000 description 8
- 239000002872 contrast media Substances 0.000 description 8
- 238000010586 diagram Methods 0.000 description 8
- 210000000056 organ Anatomy 0.000 description 6
- 210000001367 artery Anatomy 0.000 description 5
- 210000004556 brain Anatomy 0.000 description 3
- 230000000004 hemodynamic effect Effects 0.000 description 3
- 239000000284 extract Substances 0.000 description 2
- 238000012905 input function Methods 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 238000009206 nuclear medicine Methods 0.000 description 2
- 230000004044 response Effects 0.000 description 2
- 238000003325 tomography Methods 0.000 description 2
- 229910052688 Gadolinium Inorganic materials 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 239000008280 blood Substances 0.000 description 1
- 210000004369 blood Anatomy 0.000 description 1
- 210000004204 blood vessel Anatomy 0.000 description 1
- 230000000593 degrading effect Effects 0.000 description 1
- 230000002542 deteriorative effect Effects 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000010230 functional analysis Methods 0.000 description 1
- UIWYJDYFSGRHKR-UHFFFAOYSA-N gadolinium atom Chemical compound [Gd] UIWYJDYFSGRHKR-UHFFFAOYSA-N 0.000 description 1
- 230000003902 lesion Effects 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000002600 positron emission tomography Methods 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Images
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
- A61B6/50—Clinical applications
- A61B6/507—Clinical applications involving determination of haemodynamic parameters, e.g. perfusion CT
Description
【０００１】
【発明の属する技術分野】
本発明は、血流動態解析装置及び方法に係り、特にＸ線や電磁波を用いて被検体の断層像から血流動態の情報を得る血流動態解析装置及び方法に関する。
【０００２】
【従来の技術】
血流動態を計測するコンピュータ断層診断装置には、陽電子放射断層撮影装置（Positron Emission Tomography：PET ）や単光子放射ＣＴ（Single Photon Emission CT ：SPECT ）などの核医学装置がある。核医学装置による血流動態計測では、放射性核種を被検体内に注入し、核種から放射された放射線をシンチレーションカメラで計測することで、横断面についての放射性核種分布像が得られる。この放射性核種分布像を解析することにより、血流量や平均通過時間に代表される血流動態の情報を得ていた。
【０００３】
血流動態を計測する他のコンピュータ断層診断装置には、Ｘ線ＣＴ装置がある。単純Ｘ線ＣＴ像により病変の超早期相を診断することが困難な場合に、Ｘ線ＣＴ装置による血流動態計測では、コントラストを付けるために主にヨード系の造影剤を用いていた。造影剤の濃度とＣＴ値とは比例関係にあるため、造影剤を被検体に注入しダイナミックスキャンを行うことで、撮影断面の時間変化の情報が得られる。この撮影断面の時間変化から、各組織のＣＴ値の時間変化を示す時間−濃度曲線が得られ、この時間−濃度曲線を解析することにより、血流量や平均通過時間に代表される血流動態の情報を得ていた。
【０００４】
血流動態を計測する更に他のコンピュータ断層診断装置には、ＭＲＩ（Magnetic Resonance Imaging）装置がある。高磁場ＭＲＩ装置により非造影で診断する方法もあるが、高磁場ＭＲＩ装置によりコントラストの良好な画像を得る場合、低磁場ＭＲＩ装置により診断する場合においては、主にガドリニウム系の造影剤を用いていた。ＭＲＩ装置についても、撮影断面の時間変化の情報から、各組織のＭＲ値の時間変化を示す時間−濃度曲線が得られ、この時間−濃度曲線を解析することにより、血流量や平均通過時間に代表される血流動態の情報を得ていた。
【０００５】
Ｘ線ＣＴ装置やＭＲＩ装置等の画像診断装置により、脳などの器官の血流動態を調べる場合に、コントラストの良好な投影像を得るために造影剤を用いるが、この造影剤は被検体内に急速に注入する必要がある。造影剤の注入量には限度があるため、組織のＣＴ値やＭＲ信号値の変化は、ノイズを打ち消す程度に大きくならない。このため、ＳＮ比（signal-to-noise ratio : 信号と雑音との振幅の比）を上げるためにＣＴ画像やＭＲＩ画像に対して一定のマトリクスサイズでスムージング処理を行っていた。
【０００６】
高磁場ＭＲＩ装置により脳などの器官の血流動態を調べる場合、造影剤を用いずに検査することも可能である。この場合、撮影時に印加できる傾斜磁場の大きさには限度があるため、ＭＲ信号値の変化はノイズを打ち消す程度に大きくならない。このため、ＳＮ比を上げるためにＭＲＩ画像に対して一定のマトリクスサイズでスムージング処理を行っていた。
【０００７】
従来、血流動態を解析するアルゴリズムには、ガンマ関数解析法がある。この方法は、ダイナミックスキャンによって得られ、ＣＴ値やＭＲ信号値の時間変化の情報を示す時間−濃度曲線に対して、ガンマ関数でカーブフィッティングを行い、フィッティング曲線の曲線下面積やピーク値から血流量を算出していた。
【０００８】
血流動態を解析する他のアルゴリズムには、最大勾配（maximum slope ）法がある。この方法は、時間−濃度曲線の最大勾配から血流量を算出していた。最大勾配法は、ガンマ関数解析法に比べて定量性が改善される。
【０００９】
血流動態を解析する更に他のアルゴリズムには、逆たたみこみ（deconvolution ）法がある。この方法は、動脈の時間−濃度曲線と各画素の時間−濃度曲線とをdeconvolution 演算することでインパルス応答を算出し、このインパルス応答の曲線下面積やピーク値から血流量を算出していた。deconvolution 法は、ガンマ関数解析法や最大勾配法に比べて、定量性や造影レートについての問題が改善される。
【００１０】
【発明が解決しようとする課題】
ところで、機能解析を行う器官について、断層像でのコントラストの付き具合は画素毎に異なるため、時間−濃度曲線のＳＮ比も画素毎に異なる。従来、断層像のすべての画素に対して同じマトリクスサイズでスムージングを行っていた。このため、時間−濃度曲線にてＳＮ比の低い画素でのノイズを打ち消すため、スムージングのマトリクスサイズを大きくすると、ＳＮ比の高い画素での分解能が必要以上に低下するという問題があった。
【００１１】
分解能が必要以上に低下するのを防止するためには、スムージングのマトリクスサイズを小さくすればよい。しかし、マトリクスサイズを小さくした場合には、時間−濃度曲線のＳＮ比が低い組織についてノイズを打ち消すことができないという問題があった。
【００１２】
機能解析する器官の断層像についてスムージングする際、最適なマトリクスサイズは撮影条件や撮影部位等により異なる。このため、最適なマトリクスサイズを選択するには、オペレータの経験と習熟とが必要とされるという問題があった。
【００１３】
最適なマトリクスサイズが選択された場合でも、最適にスムージングを行うための重み付けは、撮影対象の性質により異なる。そのため、最適にスムージングを行うための重み付けを決定するには、オペレータの経験と習熟とが必要とされるという問題があった。
【００１４】
deconvolution 演算は、大別して行列式を用いる方法と周波数空間に変換する方法とがある。周波数空間に変換する方法を用いた場合、周波数空間でのdeconvolution 演算は、高周波数成分を必要以上に強調するため、ノイズをできる限り排除する必要がある。このため、周波数空間でdeconvolution 演算する場合には、一般に高周波数領域を抑制するフィルタ（高周波数領域制御フィルタ）によりフィルタリングを行う。フィルタリングを行う際に、最適な高周波数領域制御フィルタの形状は、撮影条件や撮影部位等に依存して決定されるため、オペレータの経験と習熟とが必要とされるという問題があった。
【００１５】
本発明はこのような事情に鑑みてなされたもので、断層像から血流動態の情報を得る場合に、機能画像の分解能を必要以上に低下させることなく、かつ良好にノイズを除去することができる血流動態解析装置及び方法を提供することを目的とする。
【００１６】
【課題を解決するための手段】
請求項１に係る血流動態解析装置は、コンピュータ断層像撮影手段で撮影した被検体の断層像を入力する断層像入力手段と、前記断層像における各画素ごとに画素値の時間変化情報を取り出す手段と、各画素ごとに取り出した時間変化情報からノイズを除去するノイズ除去手段を有し、各画素ごとにノイズ除去した時間変化情報に基づいて血流動態を解析する演算手段と、からなる血流動態解析装置において、前記ノイズ除去手段は、各画素ごとに取り出された画素値の時間変化情報と、該画素値の時間変化情報をカーブフィッティングした関数との誤差を示す値であるパラメータを算出するパラメータ算出手段と、前記算出したパラメータに対応するマトリクスサイズを有する画像フィルタを選択するフィルタ選択手段と、を有し、画素ごとに画像フィルタを選択してスムージング処理を行うことを特徴としている。
請求項２に係る血流動態解析装置は、前記フィルタ選択手段は、前記パラメータと、前記画像フィルタとを対応付けて記憶するルックアップテーブルを有し、前記算出したパラメータに基づいて前記ルックアップテーブルから対応する画像フィルタを選択することを特徴とする。
請求項３に係る血流動態解析装置は、前記パラメータは、前記画素値の時間変化情報と、該画素値の時間変化情報をカーブフィッティングした関数との誤差の絶対値の標準偏差値であることを特徴とする。
【００１７】
請求項４に係る血流動態解析装置は、コンピュータ断層像撮影手段で撮影した被検体の断層像を入力する断層像入力手段と、前記断層像における各画素ごとに画素値の時間変化情報を取り出す手段と、各画素ごとに取り出した時間変化情報の周波数空間でのノイズを除去するノイズ除去手段を有し、各画素ごとにノイズ除去した時間変化情報に基づいて血流動態を解析する演算手段と、からなる血流動態解析装置において、前記ノイズ除去手段は、各画素ごとに取り出された時間変化情報に基づいて高周波数領域制御フィルタの形状を決定するフィルタ形状決定手段を有し、画素ごとに高周波数領域制御フィルタを決定してフィルタリング処理を行うことを特徴とする。
【００１８】
請求項５に記載の血流動態解析方法は、コンピュータ断層像撮影手段で撮影した被検体の断層像の入力を受付ける工程と、前記断層像における各画素ごとに画素値の時間変化情報を取り出す工程と、各画素ごとに取り出された画素値の時間変化情報と、該画素値の時間変化情報をカーブフィッティングした関数との誤差に基づく値であるパラメータを算出する工程と、前記算出したパラメータに対応するマトリクスサイズを有する画像フィルタを選択するフィルタ選択工程と、前記フィルタ選択工程により選択された画素ごとの画像フィルタに基づいてスムージング処理を行なうことにより、各画素ごとに取り出した時間変化情報からノイズを除去する工程と、前記ノイズを除去する工程により各画素ごとにノイズ除去した時間変化情報に基づいて血流動態を解析する演算工程と、を含むことを特徴とする。
また、請求項６に係る血流動態解析方法は、コンピュータ断層像撮影手段で撮影した被検体の断層像の入力を受付ける断層像入力工程と、前記断層像における各画素ごとに画素値の時間変化情報を取り出す工程と、各画素ごとに取り出された時間変化情報に基づいて高周波数領域制御フィルタの形状を決定するフィルタ形状決定工程と、画素ごとに高周波数領域制御フィルタを決定してフィルタリング処理を行うことにより各画素ごとに取り出した時間変化情報の周波数空間でのノイズを除去する工程と、前記ノイズを除去する工程により各画素ごとにノイズ除去した時間変化情報に基づいて血流動態を解析する演算工程と、を含むことを特徴とする。
【００１９】
本発明によれば、断層像から血流動態の解析を行う場合、入力画像から血流動態の情報を解析することにより、オペレータの経験と習熟がなくても、撮影条件や撮影部位等によって異なる最適な画像フィルタ（最適なマトリクスサイズと重み付け等を有するフィルタ）を、画素毎に自動で決定することができる。
【００２０】
また、周波数空間でノイズ除去（高周波数領域制御）を行う場合、画素ごとに入力画像から血流動態の情報を解析することにより、オペレータの経験と習熟がなくても、撮影条件や撮影部位等によって異なる最適な高周波数領域制御フィルタの形状を自動で決定することができる。
【００２１】
【発明の実施の形態】
以下添付図面に従って本発明に係る血流動態解析装置及び方法の好ましい実施の形態について詳説する。
【００２２】
図１は、本発明に係る血流動態解析装置のハードウェア構成例を示すブロック図である。同図に示すように、本発明に係る血流動態解析装置は、主として断層像撮影手段で撮影した断層像を入力する断層像入力装置１と、血流動態の解析などの諸演算を行う演算装置２とから構成され、機能画像の分解能の低下を防止する。
【００２３】
断層像入力装置１と演算装置２は必ずしも別個の装置である必要はなく、一体化していてもよい。また、断層像入力装置１と演算装置２は、必ずしも断層像撮影手段と別個の装置である必要はなく、一体化していてもよい。
【００２４】
演算装置２は、断層像入力装置１から断層像の取り込みを行うインターフェイス（Ｉ／Ｆ）３と、断層像や演算結果を一時的に格納するメモリ４と、血流動態の解析などの諸演算を行う中央処理装置（ＣＰＵ）５と、演算結果や機能画像を記録するハードディスクなどの記録媒体６と、マウスやキーボードなどの外部入力装置７と、処理結果を表示する表示装置８と、各構成要素を接続する共通バス１０とから構成される。
【００２５】
図２は、本発明に係る血流動態解析装置における断層像の入力から機能画像の表示までの処理を示すフローチャートである。図２に示すように、まず、断層像入力装置１から血流解析を行う被検体の断層像を選択し、選択された断層像を演算装置２に入力する（ステップＳ２０１）。
【００２６】
入力された断層像はメモリ４に一時格納されるか、あるいは記録媒体６に保存される。また、断層像入力装置１と演算装置２とが一体化して構成される場合には、血流解析を行う被検体の断層像を外部入力装置７により選択し、選択された断層像を記録媒体６から読み出し、メモリ４に一時格納する。続いて、中央処理装置５により格納された断層像から、画素値が時間変化する情報を表す時間−濃度曲線を画素毎に取り出す（ステップＳ２０２）。次に、取り出された時間−濃度曲線から、Ｐｅａｋ値やＳＤ値（standard deviation：標準偏差）などのマトリクスサイズ決定パラメータを算出する（ステップＳ２０３）。尚、マトリクスサイズ決定パラメータの算出方法については後述する。
【００２７】
次に、中央処理装置５により、算出したマトリクスサイズ決定パラメータに対応して、スムージング処理の際のマトリクスサイズを決定する（ステップＳ２０４）。任意のマトリクスサイズ決定パラメータに対応するマトリクスサイズを決定するため、予めルックアップテーブル等の対応表を作成しておく。尚、対応表の作成方法については後述する。
【００２８】
次に、各断層像における局所領域のＳＤ値等の重み付け決定パラメータを算出する (ステップＳ２０５) 。尚、重み付け決定パラメータの算出方法については後述する。次に、中央処理装置５により、算出した重み付け決定パラメータから、最適な重み付けを決定する (ステップＳ２０６) 。尚、重み付けの決定方法については後述する。
【００２９】
次に、決定したマトリクスサイズおよび重み付けの画像フィルタにより、画素毎にスムージング処理を行う（ステップＳ２０７）。ステップＳ２０２からステップＳ２０７までの処理を断層像の全画素に対して行うことにより、スムージング処理を施した断層像が得られる。次に、中央処理装置５により、スムージング処理が施された断層像から、画素毎に時間−濃度曲線を取り出す（ステップＳ２０８）。次に、中央処理装置５により、この時間−濃度曲線に対して画素毎に機能情報を算出する（ステップＳ２０９）。この機能情報を解析する際のアルゴリズムは、ガンマ関数解析法、最大勾配法、deconvolution 法等の公知のアルゴリズムでよい。ステップＳ２０８とステップＳ２０９の処理は、スムージング処理後の断層像の全画素に対して行う。次に、断層像の全画素について解析された機能情報をマッピングすることで機能画像が得られ（ステップＳ２１０）、この機能画像は表示装置８に表示される（ステップＳ２１１）。
【００３０】
次に、時間−濃度曲線のＰｅａｋ値やＳＤ値などのマトリクスサイズ決定パラメータの算出方法について説明する。図３（ａ）、（ｂ）は、時間−濃度曲線のＳＤ値を算出する方法の説明図である。図３（ａ）に示すように、断層像の任意の画素から時間−濃度曲線を取り出し、この時間−濃度曲線の最大値をＰｅａｋ値とする。続いて、取り出した時間−濃度曲線に対して、ガンマ関数等の公知の関数でカーブフィッティングを行う（図３（ｂ））。尚、カーブフィッティングを行う際のアルゴリズムは、最小二乗法等の公知のアルゴリズムで良い。また、フィッティング関数により得られた最大値をＰｅａｋ値としても良い。
【００３１】
ここで、濃度時間曲線をTDC(t)、フィッティング関数をfit(t)とすると、両者の誤差の絶対値g(t)は次式で表される。
【００３２】
【数１】
g(t)のＳＤ値は次式で表される。
【００３３】
【数２】
ここで、ｇａｖ（ｔ）は、ｇ（ｔ）の平均値を表す。[数２]式により算出されるＳＤ値を、時間−濃度曲線のＳＤ値と定義する。
【００３４】
ここでは、マトリクスサイズ決定パラメータの代表例として、Ｐｅａｋ値とＳＤ値とを取り上げたが、マトリクスサイズ決定パラメータはこれらに限定されない。例えば［数１］式に示す誤差の絶対値g(t)又は時間−濃度曲線の分散、あるいはＰｅａｋ値とＳＤ値との比など、時間−濃度曲線が解析精度を確保するのに十分なＳＮ比を持つように、スムージング処理を行う際のマトリクスサイズを決定するパラメータならば、いかなるものでもよい。
【００３５】
次に、マトリクスサイズ決定パラメータとマトリクスサイズとの対応表を作成する方法について図４（ａ）〜（ｅ）を用いて説明する。
【００３６】
まず、図４（ａ）に示すように、血流動態を解析する器官と同等の大きさの水ファントム４０１中に任意の径の管４０２を配置したファントム４００を用意する。このファントム４００外壁の断層像の形状は、血流動態を解析する器官の断層像の形状と類似していることが望ましい。例えば、人間の脳の血流動態を解析する場合には、直径１６０mm〜１９０mm程度の円筒形のファントム等が適切である。
【００３７】
ファントム４００内部の模擬血管を示す管４０２に、任意の造影条件（濃度、量、造影レートなど）で造影剤を流しながら任意の撮影条件（管電流、管電圧、スライス厚、スキャン時間など）で撮影し、断層像を作成する。この断層像から時間−濃度曲線を取得し、この時間−濃度曲線に対して算出したＳＤ値をａとする（図４（ｂ））。ここでＳＤ値ａが閾値ＴＨよりも大きい場合、Ｎ×Ｎのマトリクスでスムージング処理を行う（図４（ｃ））。Ｎは、スムージング処理後の断層像のＳＤ値ｂが閾値ＴＨ以下になるような最小の整数である。閾値ＴＨは、十分な解析精度が得られる範囲内で任意に設定すればよい。これにより、マトリクスサイズ決定パラメータとマトリクスサイズとの関係の一例が得られる。同様にして、設定可能な複数の造影条件、撮影条件でファントムを撮影することにより、様々なマトリクスサイズ決定パラメータにおけるマトリクスサイズの関係が得られる（図４（ｄ））。これらの結果から、マトリクスサイズはマトリクスサイズ決定パラメータの何らかの関数で表すことができる。よって、任意のマトリクスサイズ決定パラメータに対応するマトリクスサイズを求めることができ、図４（ｅ）に示す対応表を作成することができる。この対応表を用いれば、画素毎に最適なマトリクスサイズを自動決定することが可能である。尚、本実施の形態では、便宜上、ファントム４００内部に一本の管４０２のみを配置したが、管の径や本数は任意でよい。
【００３８】
次に、重み付け決定パラメータの算出方法および重み付けの決定方法について図５（ａ）〜（ｅ）を用いて説明する。
【００３９】
例えば、図５（ａ）に示すように、断層像５００に対してＮ×Ｎのマトリクスでスムージング処理を行う場合において、スムージングを行う画素の座標を（i,j ）とおく。ここで、i-(N-1)/2 ≦i ≦i ，j-(N-1)/2 ≦j ≦j の領域（図５（ｂ））の標準偏差をＳＤ１、スムージングの重み付けをｗ１とし、i ≦i ≦i+(N-1)/2 ，j-(N-1)/2 ≦j ≦j の領域（図５（ｃ））の標準偏差をＳＤ２、スムージングの重み付けをｗ２とし、i-(N-1)/2 ≦i ≦i ，j ≦j ≦j+(N-1)/2 の領域（図５（ｄ））の標準偏差をＳＤ３、スムージング処理の重み付けをｗ３とし、i ≦i ≦i+(N-1)/2 ，j ≦j ≦j+(N-1)/2 の領域（図５（ｅ））の標準偏差をＳＤ４、スムージング処理の重み付けをｗ４とする。ここで必要以上に分解能を低下させることなくスムージング処理を行うためには、標準偏差の小さい領域の重み付けを大きく、標準偏差の大きい領域の重み付けを小さく設定すればよく、例えば次式のように、スムージング処理の重み付けを標準偏差の逆数の比に等しくなるように設定すればよい。
【００４０】
【数３】
［数３］式のように重み付けを算出することで、必要以上に分解能を低下させることなく、オペレータの経験や習熟に依らずにスムージング処理の重み付けを自動決定することが可能である。
【００４１】
また、重み付け方法は［数３］式に限らず、例えば最も標準偏差の小さい領域の重み付けを１、その他の領域の重み付けを０に決定する等の、必要以上に分解能が低下しないならば任意の重み付け方法を用いてもよい。
【００４２】
更に、本実施の形態では、便宜上、重み付け決定パラメータに標準偏差を用いたが、重み付け決定パラメータはこれに限らず、分散やスムージング処理を行う画素とマトリクス内のその他の画素との差の絶対値和等の、マトリクス内の画素値のばらつきを判定できるパラメータであれば任意のパラメータでよい。
【００４３】
更に、本実施の形態では、便宜上、標準偏差を算出する領域を図５（ｂ）から図５（ｅ）までの４つの領域に分割したが、分割の方向や形状、個数等の分割方法はこれに限らず、必要以上に分解能が低下しないならば任意の分割方法でよい。
【００４４】
次に、ステップＳ２０９にて機能情報を解析するためのアルゴリズムにdeconvolution 法を用いた場合に、高周波数領域制御フィルタの形状を自動決定する方法について図６（ａ）〜（ｃ）を用いて説明する。
【００４５】
前述のように、deconvolution 演算を周波数空間で行う場合、一般に、図７に示す形状の高周波数領域制御フィルタでフィルタリングすることにより、無用なノイズ強調を抑制する。例えば、対象組織への動脈入力関数、即ち流入動脈における時間−濃度曲線が、図６（ａ）に示すように実測されたとする。この時間−濃度曲線の平均通過時間（ＭＴＴ）は、次式で求められる。
【００４６】
【数４】
ここで、ｔは時間、ｎは計測点数、Ca(t) は時間−濃度曲線を表す。流入動脈における時間−濃度曲線のピークタイムＰＴ（図６（ａ））と、［数４］式から求められる平均通過時間から、図６（ｂ）に示すようなピークタイムＰＴを中心に平均通過時間の幅を持つ理想的な矩形の入力関数ideal ＿Ca(t) が求められる。このideal ＿Ca(t) のフーリエ変換をＦ｛ideal ＿Ca(t) ｝とすると、図６（ｃ）に示すようになる。
【００４７】
このＦ｛ideal ＿Ca(t) ｝に対して、例えば、周波数０から２π／ｋまでの成分のみをそのまま残し、２π／ｋから４π／ｋまでの成分を減弱し、４π／ｋよりも高い周波数成分はカットするようにフィルタの形状を決定する。この高周波数領域制御フィルタは、次式で与えられる。
【００４８】
【数５】
ここで、ωs は周波数制御開始周波数、ωe は周波数制御終了周波数であり、ideal ＿Ca(t) の幅、すなわちCa(t) の平均通過時間をｋとすれば、ωs ＝２π／ｋ，ωe ＝４π／ｋである。よって、流入動脈における時間−濃度曲線から、高周波領域制御フィルタの形状を決定することができる。
【００４９】
本実施の形態では、便宜上、ωs ＝２π／ｋ，ωe ＝４π／ｋとして説明したが、ωs とωe はこれに限らず、ωs ＝ｍπ／ｋ，ωe ＝ｎπ／ｋ（ｍ≦ｎ）として任意に設定してもよい。更に、本実施の形態では、流入動脈における時間−濃度曲線からωs とωe を算出したが、ωs とωe の算出方法はこれに限らず、図７のように最大周波数から任意に設定してもよい。
【００５０】
【発明の効果】
以上説明したように本発明に係る血流動態解析装置及び方法によれば、血流動態の情報を得る場合に、機能画像の分解能を必要以上に低下させることなく、画素毎に最適な画像フィルタを選択してスムージング処理を行なうことができる。
【００５１】
更に、周波数空間で高周波数領域の制御を行う場合に、最適な高周波数領域制御フィルタの形状を画素ごとに自動で決定してフィルタリング処理を行うことができる。
【図面の簡単な説明】
【図１】本発明に係る血流動態解析装置のハードウェア構成例を示すブロック図
【図２】本発明に係る血流動態解析装置における処理手順を示すフローチャート
【図３】時間−濃度曲線のＳＤ値を算出する方法を説明する図
【図４】マトリクスサイズ決定パラメータとマトリクスサイズとの対応表を作成する方法を説明する図
【図５】重み付け決定パラメータと重み付けの決定方法を説明する図
【図６】高周波数領域制御フィルタの作成方法を説明する図
【図７】高周波数領域制御フィルタの例を示す図
【符号の説明】
１…断層像入力装置、２…演算装置、３…インターフェイス（Ｉ／Ｆ）、４…メモリ、５…中央処理装置（ＣＰＵ）、６…記録媒体（ハードディスク）、７…外部入力装置、８…表示装置
【発明の属する技術分野】
本発明は、血流動態解析装置及び方法に係り、特にＸ線や電磁波を用いて被検体の断層像から血流動態の情報を得る血流動態解析装置及び方法に関する。
【０００２】
【従来の技術】
血流動態を計測するコンピュータ断層診断装置には、陽電子放射断層撮影装置（Positron Emission Tomography：PET ）や単光子放射ＣＴ（Single Photon Emission CT ：SPECT ）などの核医学装置がある。核医学装置による血流動態計測では、放射性核種を被検体内に注入し、核種から放射された放射線をシンチレーションカメラで計測することで、横断面についての放射性核種分布像が得られる。この放射性核種分布像を解析することにより、血流量や平均通過時間に代表される血流動態の情報を得ていた。
【０００３】
血流動態を計測する他のコンピュータ断層診断装置には、Ｘ線ＣＴ装置がある。単純Ｘ線ＣＴ像により病変の超早期相を診断することが困難な場合に、Ｘ線ＣＴ装置による血流動態計測では、コントラストを付けるために主にヨード系の造影剤を用いていた。造影剤の濃度とＣＴ値とは比例関係にあるため、造影剤を被検体に注入しダイナミックスキャンを行うことで、撮影断面の時間変化の情報が得られる。この撮影断面の時間変化から、各組織のＣＴ値の時間変化を示す時間−濃度曲線が得られ、この時間−濃度曲線を解析することにより、血流量や平均通過時間に代表される血流動態の情報を得ていた。
【０００４】
血流動態を計測する更に他のコンピュータ断層診断装置には、ＭＲＩ（Magnetic Resonance Imaging）装置がある。高磁場ＭＲＩ装置により非造影で診断する方法もあるが、高磁場ＭＲＩ装置によりコントラストの良好な画像を得る場合、低磁場ＭＲＩ装置により診断する場合においては、主にガドリニウム系の造影剤を用いていた。ＭＲＩ装置についても、撮影断面の時間変化の情報から、各組織のＭＲ値の時間変化を示す時間−濃度曲線が得られ、この時間−濃度曲線を解析することにより、血流量や平均通過時間に代表される血流動態の情報を得ていた。
【０００５】
Ｘ線ＣＴ装置やＭＲＩ装置等の画像診断装置により、脳などの器官の血流動態を調べる場合に、コントラストの良好な投影像を得るために造影剤を用いるが、この造影剤は被検体内に急速に注入する必要がある。造影剤の注入量には限度があるため、組織のＣＴ値やＭＲ信号値の変化は、ノイズを打ち消す程度に大きくならない。このため、ＳＮ比（signal-to-noise ratio : 信号と雑音との振幅の比）を上げるためにＣＴ画像やＭＲＩ画像に対して一定のマトリクスサイズでスムージング処理を行っていた。
【０００６】
高磁場ＭＲＩ装置により脳などの器官の血流動態を調べる場合、造影剤を用いずに検査することも可能である。この場合、撮影時に印加できる傾斜磁場の大きさには限度があるため、ＭＲ信号値の変化はノイズを打ち消す程度に大きくならない。このため、ＳＮ比を上げるためにＭＲＩ画像に対して一定のマトリクスサイズでスムージング処理を行っていた。
【０００７】
従来、血流動態を解析するアルゴリズムには、ガンマ関数解析法がある。この方法は、ダイナミックスキャンによって得られ、ＣＴ値やＭＲ信号値の時間変化の情報を示す時間−濃度曲線に対して、ガンマ関数でカーブフィッティングを行い、フィッティング曲線の曲線下面積やピーク値から血流量を算出していた。
【０００８】
血流動態を解析する他のアルゴリズムには、最大勾配（maximum slope ）法がある。この方法は、時間−濃度曲線の最大勾配から血流量を算出していた。最大勾配法は、ガンマ関数解析法に比べて定量性が改善される。
【０００９】
血流動態を解析する更に他のアルゴリズムには、逆たたみこみ（deconvolution ）法がある。この方法は、動脈の時間−濃度曲線と各画素の時間−濃度曲線とをdeconvolution 演算することでインパルス応答を算出し、このインパルス応答の曲線下面積やピーク値から血流量を算出していた。deconvolution 法は、ガンマ関数解析法や最大勾配法に比べて、定量性や造影レートについての問題が改善される。
【００１０】
【発明が解決しようとする課題】
ところで、機能解析を行う器官について、断層像でのコントラストの付き具合は画素毎に異なるため、時間−濃度曲線のＳＮ比も画素毎に異なる。従来、断層像のすべての画素に対して同じマトリクスサイズでスムージングを行っていた。このため、時間−濃度曲線にてＳＮ比の低い画素でのノイズを打ち消すため、スムージングのマトリクスサイズを大きくすると、ＳＮ比の高い画素での分解能が必要以上に低下するという問題があった。
【００１１】
分解能が必要以上に低下するのを防止するためには、スムージングのマトリクスサイズを小さくすればよい。しかし、マトリクスサイズを小さくした場合には、時間−濃度曲線のＳＮ比が低い組織についてノイズを打ち消すことができないという問題があった。
【００１２】
機能解析する器官の断層像についてスムージングする際、最適なマトリクスサイズは撮影条件や撮影部位等により異なる。このため、最適なマトリクスサイズを選択するには、オペレータの経験と習熟とが必要とされるという問題があった。
【００１３】
最適なマトリクスサイズが選択された場合でも、最適にスムージングを行うための重み付けは、撮影対象の性質により異なる。そのため、最適にスムージングを行うための重み付けを決定するには、オペレータの経験と習熟とが必要とされるという問題があった。
【００１４】
deconvolution 演算は、大別して行列式を用いる方法と周波数空間に変換する方法とがある。周波数空間に変換する方法を用いた場合、周波数空間でのdeconvolution 演算は、高周波数成分を必要以上に強調するため、ノイズをできる限り排除する必要がある。このため、周波数空間でdeconvolution 演算する場合には、一般に高周波数領域を抑制するフィルタ（高周波数領域制御フィルタ）によりフィルタリングを行う。フィルタリングを行う際に、最適な高周波数領域制御フィルタの形状は、撮影条件や撮影部位等に依存して決定されるため、オペレータの経験と習熟とが必要とされるという問題があった。
【００１５】
本発明はこのような事情に鑑みてなされたもので、断層像から血流動態の情報を得る場合に、機能画像の分解能を必要以上に低下させることなく、かつ良好にノイズを除去することができる血流動態解析装置及び方法を提供することを目的とする。
【００１６】
【課題を解決するための手段】
請求項１に係る血流動態解析装置は、コンピュータ断層像撮影手段で撮影した被検体の断層像を入力する断層像入力手段と、前記断層像における各画素ごとに画素値の時間変化情報を取り出す手段と、各画素ごとに取り出した時間変化情報からノイズを除去するノイズ除去手段を有し、各画素ごとにノイズ除去した時間変化情報に基づいて血流動態を解析する演算手段と、からなる血流動態解析装置において、前記ノイズ除去手段は、各画素ごとに取り出された画素値の時間変化情報と、該画素値の時間変化情報をカーブフィッティングした関数との誤差を示す値であるパラメータを算出するパラメータ算出手段と、前記算出したパラメータに対応するマトリクスサイズを有する画像フィルタを選択するフィルタ選択手段と、を有し、画素ごとに画像フィルタを選択してスムージング処理を行うことを特徴としている。
請求項２に係る血流動態解析装置は、前記フィルタ選択手段は、前記パラメータと、前記画像フィルタとを対応付けて記憶するルックアップテーブルを有し、前記算出したパラメータに基づいて前記ルックアップテーブルから対応する画像フィルタを選択することを特徴とする。
請求項３に係る血流動態解析装置は、前記パラメータは、前記画素値の時間変化情報と、該画素値の時間変化情報をカーブフィッティングした関数との誤差の絶対値の標準偏差値であることを特徴とする。
【００１７】
請求項４に係る血流動態解析装置は、コンピュータ断層像撮影手段で撮影した被検体の断層像を入力する断層像入力手段と、前記断層像における各画素ごとに画素値の時間変化情報を取り出す手段と、各画素ごとに取り出した時間変化情報の周波数空間でのノイズを除去するノイズ除去手段を有し、各画素ごとにノイズ除去した時間変化情報に基づいて血流動態を解析する演算手段と、からなる血流動態解析装置において、前記ノイズ除去手段は、各画素ごとに取り出された時間変化情報に基づいて高周波数領域制御フィルタの形状を決定するフィルタ形状決定手段を有し、画素ごとに高周波数領域制御フィルタを決定してフィルタリング処理を行うことを特徴とする。
【００１８】
請求項５に記載の血流動態解析方法は、コンピュータ断層像撮影手段で撮影した被検体の断層像の入力を受付ける工程と、前記断層像における各画素ごとに画素値の時間変化情報を取り出す工程と、各画素ごとに取り出された画素値の時間変化情報と、該画素値の時間変化情報をカーブフィッティングした関数との誤差に基づく値であるパラメータを算出する工程と、前記算出したパラメータに対応するマトリクスサイズを有する画像フィルタを選択するフィルタ選択工程と、前記フィルタ選択工程により選択された画素ごとの画像フィルタに基づいてスムージング処理を行なうことにより、各画素ごとに取り出した時間変化情報からノイズを除去する工程と、前記ノイズを除去する工程により各画素ごとにノイズ除去した時間変化情報に基づいて血流動態を解析する演算工程と、を含むことを特徴とする。
また、請求項６に係る血流動態解析方法は、コンピュータ断層像撮影手段で撮影した被検体の断層像の入力を受付ける断層像入力工程と、前記断層像における各画素ごとに画素値の時間変化情報を取り出す工程と、各画素ごとに取り出された時間変化情報に基づいて高周波数領域制御フィルタの形状を決定するフィルタ形状決定工程と、画素ごとに高周波数領域制御フィルタを決定してフィルタリング処理を行うことにより各画素ごとに取り出した時間変化情報の周波数空間でのノイズを除去する工程と、前記ノイズを除去する工程により各画素ごとにノイズ除去した時間変化情報に基づいて血流動態を解析する演算工程と、を含むことを特徴とする。
【００１９】
本発明によれば、断層像から血流動態の解析を行う場合、入力画像から血流動態の情報を解析することにより、オペレータの経験と習熟がなくても、撮影条件や撮影部位等によって異なる最適な画像フィルタ（最適なマトリクスサイズと重み付け等を有するフィルタ）を、画素毎に自動で決定することができる。
【００２０】
また、周波数空間でノイズ除去（高周波数領域制御）を行う場合、画素ごとに入力画像から血流動態の情報を解析することにより、オペレータの経験と習熟がなくても、撮影条件や撮影部位等によって異なる最適な高周波数領域制御フィルタの形状を自動で決定することができる。
【００２１】
【発明の実施の形態】
以下添付図面に従って本発明に係る血流動態解析装置及び方法の好ましい実施の形態について詳説する。
【００２２】
図１は、本発明に係る血流動態解析装置のハードウェア構成例を示すブロック図である。同図に示すように、本発明に係る血流動態解析装置は、主として断層像撮影手段で撮影した断層像を入力する断層像入力装置１と、血流動態の解析などの諸演算を行う演算装置２とから構成され、機能画像の分解能の低下を防止する。
【００２３】
断層像入力装置１と演算装置２は必ずしも別個の装置である必要はなく、一体化していてもよい。また、断層像入力装置１と演算装置２は、必ずしも断層像撮影手段と別個の装置である必要はなく、一体化していてもよい。
【００２４】
演算装置２は、断層像入力装置１から断層像の取り込みを行うインターフェイス（Ｉ／Ｆ）３と、断層像や演算結果を一時的に格納するメモリ４と、血流動態の解析などの諸演算を行う中央処理装置（ＣＰＵ）５と、演算結果や機能画像を記録するハードディスクなどの記録媒体６と、マウスやキーボードなどの外部入力装置７と、処理結果を表示する表示装置８と、各構成要素を接続する共通バス１０とから構成される。
【００２５】
図２は、本発明に係る血流動態解析装置における断層像の入力から機能画像の表示までの処理を示すフローチャートである。図２に示すように、まず、断層像入力装置１から血流解析を行う被検体の断層像を選択し、選択された断層像を演算装置２に入力する（ステップＳ２０１）。
【００２６】
入力された断層像はメモリ４に一時格納されるか、あるいは記録媒体６に保存される。また、断層像入力装置１と演算装置２とが一体化して構成される場合には、血流解析を行う被検体の断層像を外部入力装置７により選択し、選択された断層像を記録媒体６から読み出し、メモリ４に一時格納する。続いて、中央処理装置５により格納された断層像から、画素値が時間変化する情報を表す時間−濃度曲線を画素毎に取り出す（ステップＳ２０２）。次に、取り出された時間−濃度曲線から、Ｐｅａｋ値やＳＤ値（standard deviation：標準偏差）などのマトリクスサイズ決定パラメータを算出する（ステップＳ２０３）。尚、マトリクスサイズ決定パラメータの算出方法については後述する。
【００２７】
次に、中央処理装置５により、算出したマトリクスサイズ決定パラメータに対応して、スムージング処理の際のマトリクスサイズを決定する（ステップＳ２０４）。任意のマトリクスサイズ決定パラメータに対応するマトリクスサイズを決定するため、予めルックアップテーブル等の対応表を作成しておく。尚、対応表の作成方法については後述する。
【００２８】
次に、各断層像における局所領域のＳＤ値等の重み付け決定パラメータを算出する (ステップＳ２０５) 。尚、重み付け決定パラメータの算出方法については後述する。次に、中央処理装置５により、算出した重み付け決定パラメータから、最適な重み付けを決定する (ステップＳ２０６) 。尚、重み付けの決定方法については後述する。
【００２９】
次に、決定したマトリクスサイズおよび重み付けの画像フィルタにより、画素毎にスムージング処理を行う（ステップＳ２０７）。ステップＳ２０２からステップＳ２０７までの処理を断層像の全画素に対して行うことにより、スムージング処理を施した断層像が得られる。次に、中央処理装置５により、スムージング処理が施された断層像から、画素毎に時間−濃度曲線を取り出す（ステップＳ２０８）。次に、中央処理装置５により、この時間−濃度曲線に対して画素毎に機能情報を算出する（ステップＳ２０９）。この機能情報を解析する際のアルゴリズムは、ガンマ関数解析法、最大勾配法、deconvolution 法等の公知のアルゴリズムでよい。ステップＳ２０８とステップＳ２０９の処理は、スムージング処理後の断層像の全画素に対して行う。次に、断層像の全画素について解析された機能情報をマッピングすることで機能画像が得られ（ステップＳ２１０）、この機能画像は表示装置８に表示される（ステップＳ２１１）。
【００３０】
次に、時間−濃度曲線のＰｅａｋ値やＳＤ値などのマトリクスサイズ決定パラメータの算出方法について説明する。図３（ａ）、（ｂ）は、時間−濃度曲線のＳＤ値を算出する方法の説明図である。図３（ａ）に示すように、断層像の任意の画素から時間−濃度曲線を取り出し、この時間−濃度曲線の最大値をＰｅａｋ値とする。続いて、取り出した時間−濃度曲線に対して、ガンマ関数等の公知の関数でカーブフィッティングを行う（図３（ｂ））。尚、カーブフィッティングを行う際のアルゴリズムは、最小二乗法等の公知のアルゴリズムで良い。また、フィッティング関数により得られた最大値をＰｅａｋ値としても良い。
【００３１】
ここで、濃度時間曲線をTDC(t)、フィッティング関数をfit(t)とすると、両者の誤差の絶対値g(t)は次式で表される。
【００３２】
【数１】
【００３３】
【数２】
【００３４】
ここでは、マトリクスサイズ決定パラメータの代表例として、Ｐｅａｋ値とＳＤ値とを取り上げたが、マトリクスサイズ決定パラメータはこれらに限定されない。例えば［数１］式に示す誤差の絶対値g(t)又は時間−濃度曲線の分散、あるいはＰｅａｋ値とＳＤ値との比など、時間−濃度曲線が解析精度を確保するのに十分なＳＮ比を持つように、スムージング処理を行う際のマトリクスサイズを決定するパラメータならば、いかなるものでもよい。
【００３５】
次に、マトリクスサイズ決定パラメータとマトリクスサイズとの対応表を作成する方法について図４（ａ）〜（ｅ）を用いて説明する。
【００３６】
まず、図４（ａ）に示すように、血流動態を解析する器官と同等の大きさの水ファントム４０１中に任意の径の管４０２を配置したファントム４００を用意する。このファントム４００外壁の断層像の形状は、血流動態を解析する器官の断層像の形状と類似していることが望ましい。例えば、人間の脳の血流動態を解析する場合には、直径１６０mm〜１９０mm程度の円筒形のファントム等が適切である。
【００３７】
ファントム４００内部の模擬血管を示す管４０２に、任意の造影条件（濃度、量、造影レートなど）で造影剤を流しながら任意の撮影条件（管電流、管電圧、スライス厚、スキャン時間など）で撮影し、断層像を作成する。この断層像から時間−濃度曲線を取得し、この時間−濃度曲線に対して算出したＳＤ値をａとする（図４（ｂ））。ここでＳＤ値ａが閾値ＴＨよりも大きい場合、Ｎ×Ｎのマトリクスでスムージング処理を行う（図４（ｃ））。Ｎは、スムージング処理後の断層像のＳＤ値ｂが閾値ＴＨ以下になるような最小の整数である。閾値ＴＨは、十分な解析精度が得られる範囲内で任意に設定すればよい。これにより、マトリクスサイズ決定パラメータとマトリクスサイズとの関係の一例が得られる。同様にして、設定可能な複数の造影条件、撮影条件でファントムを撮影することにより、様々なマトリクスサイズ決定パラメータにおけるマトリクスサイズの関係が得られる（図４（ｄ））。これらの結果から、マトリクスサイズはマトリクスサイズ決定パラメータの何らかの関数で表すことができる。よって、任意のマトリクスサイズ決定パラメータに対応するマトリクスサイズを求めることができ、図４（ｅ）に示す対応表を作成することができる。この対応表を用いれば、画素毎に最適なマトリクスサイズを自動決定することが可能である。尚、本実施の形態では、便宜上、ファントム４００内部に一本の管４０２のみを配置したが、管の径や本数は任意でよい。
【００３８】
次に、重み付け決定パラメータの算出方法および重み付けの決定方法について図５（ａ）〜（ｅ）を用いて説明する。
【００３９】
例えば、図５（ａ）に示すように、断層像５００に対してＮ×Ｎのマトリクスでスムージング処理を行う場合において、スムージングを行う画素の座標を（i,j ）とおく。ここで、i-(N-1)/2 ≦i ≦i ，j-(N-1)/2 ≦j ≦j の領域（図５（ｂ））の標準偏差をＳＤ１、スムージングの重み付けをｗ１とし、i ≦i ≦i+(N-1)/2 ，j-(N-1)/2 ≦j ≦j の領域（図５（ｃ））の標準偏差をＳＤ２、スムージングの重み付けをｗ２とし、i-(N-1)/2 ≦i ≦i ，j ≦j ≦j+(N-1)/2 の領域（図５（ｄ））の標準偏差をＳＤ３、スムージング処理の重み付けをｗ３とし、i ≦i ≦i+(N-1)/2 ，j ≦j ≦j+(N-1)/2 の領域（図５（ｅ））の標準偏差をＳＤ４、スムージング処理の重み付けをｗ４とする。ここで必要以上に分解能を低下させることなくスムージング処理を行うためには、標準偏差の小さい領域の重み付けを大きく、標準偏差の大きい領域の重み付けを小さく設定すればよく、例えば次式のように、スムージング処理の重み付けを標準偏差の逆数の比に等しくなるように設定すればよい。
【００４０】
【数３】
【００４１】
また、重み付け方法は［数３］式に限らず、例えば最も標準偏差の小さい領域の重み付けを１、その他の領域の重み付けを０に決定する等の、必要以上に分解能が低下しないならば任意の重み付け方法を用いてもよい。
【００４２】
更に、本実施の形態では、便宜上、重み付け決定パラメータに標準偏差を用いたが、重み付け決定パラメータはこれに限らず、分散やスムージング処理を行う画素とマトリクス内のその他の画素との差の絶対値和等の、マトリクス内の画素値のばらつきを判定できるパラメータであれば任意のパラメータでよい。
【００４３】
更に、本実施の形態では、便宜上、標準偏差を算出する領域を図５（ｂ）から図５（ｅ）までの４つの領域に分割したが、分割の方向や形状、個数等の分割方法はこれに限らず、必要以上に分解能が低下しないならば任意の分割方法でよい。
【００４４】
次に、ステップＳ２０９にて機能情報を解析するためのアルゴリズムにdeconvolution 法を用いた場合に、高周波数領域制御フィルタの形状を自動決定する方法について図６（ａ）〜（ｃ）を用いて説明する。
【００４５】
前述のように、deconvolution 演算を周波数空間で行う場合、一般に、図７に示す形状の高周波数領域制御フィルタでフィルタリングすることにより、無用なノイズ強調を抑制する。例えば、対象組織への動脈入力関数、即ち流入動脈における時間−濃度曲線が、図６（ａ）に示すように実測されたとする。この時間−濃度曲線の平均通過時間（ＭＴＴ）は、次式で求められる。
【００４６】
【数４】
【００４７】
このＦ｛ideal ＿Ca(t) ｝に対して、例えば、周波数０から２π／ｋまでの成分のみをそのまま残し、２π／ｋから４π／ｋまでの成分を減弱し、４π／ｋよりも高い周波数成分はカットするようにフィルタの形状を決定する。この高周波数領域制御フィルタは、次式で与えられる。
【００４８】
【数５】
【００４９】
本実施の形態では、便宜上、ωs ＝２π／ｋ，ωe ＝４π／ｋとして説明したが、ωs とωe はこれに限らず、ωs ＝ｍπ／ｋ，ωe ＝ｎπ／ｋ（ｍ≦ｎ）として任意に設定してもよい。更に、本実施の形態では、流入動脈における時間−濃度曲線からωs とωe を算出したが、ωs とωe の算出方法はこれに限らず、図７のように最大周波数から任意に設定してもよい。
【００５０】
【発明の効果】
以上説明したように本発明に係る血流動態解析装置及び方法によれば、血流動態の情報を得る場合に、機能画像の分解能を必要以上に低下させることなく、画素毎に最適な画像フィルタを選択してスムージング処理を行なうことができる。
【００５１】
更に、周波数空間で高周波数領域の制御を行う場合に、最適な高周波数領域制御フィルタの形状を画素ごとに自動で決定してフィルタリング処理を行うことができる。
【図面の簡単な説明】
【図１】本発明に係る血流動態解析装置のハードウェア構成例を示すブロック図
【図２】本発明に係る血流動態解析装置における処理手順を示すフローチャート
【図３】時間−濃度曲線のＳＤ値を算出する方法を説明する図
【図４】マトリクスサイズ決定パラメータとマトリクスサイズとの対応表を作成する方法を説明する図
【図５】重み付け決定パラメータと重み付けの決定方法を説明する図
【図６】高周波数領域制御フィルタの作成方法を説明する図
【図７】高周波数領域制御フィルタの例を示す図
【符号の説明】
１…断層像入力装置、２…演算装置、３…インターフェイス（Ｉ／Ｆ）、４…メモリ、５…中央処理装置（ＣＰＵ）、６…記録媒体（ハードディスク）、７…外部入力装置、８…表示装置
Claims (6)
- コンピュータ断層像撮影手段で撮影した被検体の断層像を入力する断層像入力手段と、前記断層像における各画素ごとに画素値の時間変化情報を取り出す手段と、各画素ごとに取り出した時間変化情報からノイズを除去するノイズ除去手段を有し、各画素ごとにノイズ除去した時間変化情報に基づいて血流動態を解析する演算手段と、からなる血流動態解析装置において、
前記ノイズ除去手段は、各画素ごとに取り出された画素値の時間変化情報と、該画素値の時間変化情報をカーブフィッティングした関数との誤差に基づく値であるパラメータを算出するパラメータ算出手段と、前記算出したパラメータに対応するマトリクスサイズを有する画像フィルタを選択するフィルタ選択手段と、を有し、画素ごとに画像フィルタを選択してスムージング処理を行うことを特徴とする血流動態解析装置。 - 前記フィルタ選択手段は、前記パラメータと、前記画像フィルタとを対応付けて記憶するルックアップテーブルを有し、前記算出したパラメータに基づいて前記ルックアップテーブルから対応する画像フィルタを選択することを特徴とする請求項１の血流動態解析装置。
- 前記パラメータは、前記画素値の時間変化情報と、該画素値の時間変化情報をカーブフィッティングした関数との誤差の絶対値の標準偏差値であることを特徴とする請求項１又は２に記載の血流動態解析装置。
- コンピュータ断層像撮影手段で撮影した被検体の断層像を入力する断層像入力手段と、前記断層像における各画素ごとに画素値の時間変化情報を取り出す手段と、各画素ごとに取り出した時間変化情報の周波数空間でのノイズを除去するノイズ除去手段を有し、各画素ごとにノイズ除去した時間変化情報に基づいて血流動態を解析する演算手段と、からなる血流動態解析装置において、
前記ノイズ除去手段は、各画素ごとに取り出された時間変化情報に基づいて高周波数領域制御フィルタの形状を決定するフィルタ形状決定手段を有し、画素ごとに高周波数領域制御フィルタを決定してフィルタリング処理を行うことを特徴とする血流動態解析装置。 - コンピュータ断層像撮影手段で撮影した被検体の断層像の入力を受付ける工程と、
前記断層像における各画素ごとに画素値の時間変化情報を取り出す工程と、
各画素ごとに取り出された画素値の時間変化情報と、該画素値の時間変化情報をカーブフィッティングした関数との誤差に基づく値であるパラメータを算出する工程と、
前記算出したパラメータに対応するマトリクスサイズを有する画像フィルタを選択するフィルタ選択工程と、
前記フィルタ選択工程により選択された画素ごとの画像フィルタに基づいてスムージング処理を行なうことにより、各画素ごとに取り出した時間変化情報からノイズを除去する工程と、
前記ノイズを除去する工程により各画素ごとにノイズ除去した時間変化情報に基づいて血流動態を解析する演算工程と、
を含むことを特徴とする血流動態解析方法。 - コンピュータ断層像撮影手段で撮影した被検体の断層像の入力を受付ける断層像入力工程と、
前記断層像における各画素ごとに画素値の時間変化情報を取り出す工程と、
各画素ごとに取り出された時間変化情報に基づいて高周波数領域制御フィルタの形状を決定するフィルタ形状決定工程と、
画素ごとに高周波数領域制御フィルタを決定してフィルタリング処理を行うことにより各画素ごとに取り出した時間変化情報の周波数空間でのノイズを除去する工程と、
前記ノイズを除去する工程により各画素ごとにノイズ除去した時間変化情報に基づいて血流動態を解析する演算工程と、
を含むことを特徴とする血流動態解析方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2002025309A JP4193168B2 (ja) | 2002-02-01 | 2002-02-01 | 血流動態解析装置及び方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2002025309A JP4193168B2 (ja) | 2002-02-01 | 2002-02-01 | 血流動態解析装置及び方法 |
Publications (3)
| Publication Number | Publication Date |
|---|---|
| JP2003225234A JP2003225234A (ja) | 2003-08-12 |
| JP2003225234A5 JP2003225234A5 (ja) | 2005-08-04 |
| JP4193168B2 true JP4193168B2 (ja) | 2008-12-10 |
Family
ID=27747497
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2002025309A Expired - Lifetime JP4193168B2 (ja) | 2002-02-01 | 2002-02-01 | 血流動態解析装置及び方法 |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP4193168B2 (ja) |
Families Citing this family (26)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| GB0326381D0 (en) * | 2003-11-12 | 2003-12-17 | Inst Of Cancer Res The | A method and means for image processing |
| EP1585058B1 (en) | 2004-03-31 | 2008-12-10 | Kabushiki Kaisha Toshiba | Medical image processing apparatus and method of processing medical image |
| US8295914B2 (en) * | 2004-11-16 | 2012-10-23 | Medrad, Inc. | Systems and methods of determining patient transfer functions and modeling patient response to a pharmaceutical injection |
| DK2990073T3 (en) | 2004-11-24 | 2018-08-13 | Bayer Healthcare Llc | DEVICES AND SYSTEMS FOR DELIVERING FLUIDS |
| JP4804039B2 (ja) * | 2005-05-27 | 2011-10-26 | 株式会社日立メディコ | 血流動態解析装置、ｘ線ｃｔ装置、ｍｒｉ装置、及び血流動態解析プログラム |
| JP4707471B2 (ja) * | 2005-06-10 | 2011-06-22 | 株式会社日立メディコ | 画像処理プログラム、装置及び方法 |
| CN101351155B (zh) * | 2006-01-05 | 2013-02-20 | 皇家飞利浦电子股份有限公司 | 可调体模 |
| GB2459075B (en) * | 2006-06-02 | 2010-12-15 | Siemens Molecular Imaging Ltd | Estimation of blood input function for functional medical scans |
| EP2097835B1 (en) | 2006-12-29 | 2018-05-30 | Bayer Healthcare LLC | Patient-based parameter generation systems for medical injection procedures |
| WO2009145076A1 (ja) * | 2008-05-28 | 2009-12-03 | 株式会社 日立メディコ | 画像処理装置、画像処理方法、及び画像処理プログラム |
| US9421330B2 (en) | 2008-11-03 | 2016-08-23 | Bayer Healthcare Llc | Mitigation of contrast-induced nephropathy |
| JP4696278B2 (ja) * | 2009-09-14 | 2011-06-08 | 国立大学法人 東京大学 | 領域分割画像生成方法、領域分割画像生成装置及びコンピュータプログラム |
| US9968309B2 (en) | 2009-12-08 | 2018-05-15 | Koninklijke Philips N.V. | Method and a correction system for correcting tracer-uptake measurements |
| JP5566174B2 (ja) * | 2010-04-27 | 2014-08-06 | ジーイー・メディカル・システムズ・グローバル・テクノロジー・カンパニー・エルエルシー | 磁気共鳴イメージング装置およびプログラム |
| RU2434288C1 (ru) * | 2010-06-08 | 2011-11-20 | Закрытое Акционерное Общество "Импульс" | Способ коррекции цифровых изображений |
| EP2585116A4 (en) | 2010-06-24 | 2017-03-29 | Bayer Healthcare LLC | Modeling of pharmaceutical propagation and parameter generation for injection protocols |
| JP6081139B2 (ja) * | 2011-12-26 | 2017-02-15 | 東芝メディカルシステムズ株式会社 | 超音波診断装置、医用画像処理装置、及び医用画像処理方法 |
| CN111528872A (zh) | 2012-05-14 | 2020-08-14 | 拜耳医药保健有限公司 | 用于基于x射线管电压确定药用流体注射协议的系统和方法 |
| US9555379B2 (en) | 2013-03-13 | 2017-01-31 | Bayer Healthcare Llc | Fluid path set with turbulent mixing chamber, backflow compensator |
| WO2017152036A1 (en) | 2016-03-03 | 2017-09-08 | Bayer Healthcare Llc | System and method for improved fluid delivery in multi-fluid injector systems |
| AU2018323442A1 (en) | 2017-08-31 | 2020-01-23 | Bayer Healthcare Llc | Fluid path impedance assessment for improving fluid delivery performance |
| US11598664B2 (en) | 2017-08-31 | 2023-03-07 | Bayer Healthcare Llc | Injector pressure calibration system and method |
| US11786652B2 (en) | 2017-08-31 | 2023-10-17 | Bayer Healthcare Llc | System and method for drive member position and fluid injector system mechanical calibration |
| CA3068544A1 (en) | 2017-08-31 | 2019-03-07 | Bayer Healthcare Llc | Method for dynamic pressure control in a fluid injector system |
| WO2019046267A1 (en) | 2017-08-31 | 2019-03-07 | Bayer Healthcare Llc | SYSTEM AND METHOD FOR VOLUME COMPENSATION OF FLUID INJECTOR SYSTEM |
| JP7355525B2 (ja) * | 2019-05-23 | 2023-10-03 | 東芝Ｉｔコントロールシステム株式会社 | 放射線検査装置 |
-
2002
- 2002-02-01 JP JP2002025309A patent/JP4193168B2/ja not_active Expired - Lifetime
Also Published As
| Publication number | Publication date |
|---|---|
| JP2003225234A (ja) | 2003-08-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP4193168B2 (ja) | 血流動態解析装置及び方法 | |
| JP6559678B2 (ja) | スペクトル画像データ処理 | |
| JP5893623B2 (ja) | データ・セットにおける異常検出方法およびシステム | |
| JP5081390B2 (ja) | 腫瘍量を監視する方法及びシステム | |
| US20050113679A1 (en) | Method and apparatus for segmenting structure in CT angiography | |
| JP4013128B2 (ja) | 血流動態解析装置、方法、及び画像診断装置 | |
| WO2010055405A2 (en) | Method and system for mapping tissue status of acute stroke | |
| JP5658686B2 (ja) | 貫壁性のかん流の勾配の画像分析 | |
| WO2005046478A1 (ja) | 画像処理方法、画像処理装置、医用画像診断支援システム、及び時間軸方向フィルタリング方法 | |
| JP4804039B2 (ja) | 血流動態解析装置、ｘ線ｃｔ装置、ｍｒｉ装置、及び血流動態解析プログラム | |
| JP2008529639A (ja) | 検査装置、画像処理デバイス、注目対象物を検査装置で検査する方法、コンピュータ可読媒体及びプログラム要素 | |
| EP3455824B1 (en) | Noise reduction in image data | |
| JP4589114B2 (ja) | 画像データ情報の表示 | |
| Chen et al. | Characterization of the image-derived carotid artery input function using independent component analysis for the quantitation of [18F] fluorodeoxyglucose positron emission tomography images | |
| Kang et al. | An anatomic coordinate system of the femoral neck for highly reproducible BMD measurements using 3D QCT | |
| Pahn et al. | Toward standardized quantitative image quality (IQ) assessment in computed tomography (CT): A comprehensive framework for automated and comparative IQ analysis based on ICRU Report 87 | |
| JP2018538032A (ja) | 機能的医用撮像における定量的マップ生成の自動最適化方法 | |
| JP2006500099A5 (ja) | ||
| JP6564075B2 (ja) | 医用画像を表示するための伝達関数の選択 | |
| US10043268B2 (en) | Medical image processing apparatus and method to generate and display third parameters based on first and second images | |
| JP5042533B2 (ja) | 医用画像表示装置 | |
| JP2003299646A (ja) | 画像解析装置 | |
| Elfarnawany et al. | Improved objective selection of power Doppler wall-filter cut-off velocity for accurate vascular quantification | |
| JP6898047B2 (ja) | 時変データの定量的評価 | |
| US11704795B2 (en) | Quality-driven image processing |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20050104 |
|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20050104 |
|
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20080401 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20080415 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20080611 |
|
| 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: 20080829 |
|
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20080911 |
|
| R150 | Certificate of patent or registration of utility model |
Ref document number: 4193168 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20111003 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20111003 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20121003 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20121003 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20131003 Year of fee payment: 5 |
|
| 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 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| 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 |