JP4363833B2

JP4363833B2 - 局所血流動態に関するインデックスを演算する方法及び装置

Info

Publication number: JP4363833B2
Application number: JP2002284450A
Authority: JP
Inventors: 恭二郎南部; 佳弘池田
Original assignee: Toshiba Corp
Current assignee: Toshiba Corp
Priority date: 2001-10-16
Filing date: 2002-09-27
Publication date: 2009-11-11
Anticipated expiration: 2022-09-27
Also published as: JP2003190148A

Description

【０００１】
【発明の属する技術分野】
本発明は、脳組織等の局所血流動態に関するインデックスの演算方法及び演算装置に関する。
【０００２】
【従来の技術】
Ｘ線ＣＴ検査では、単純ＣＴ像から形態情報を、また造影ＣＴによるダイナミックスキャンで病巣の周りの血流の動態情報をそれぞれ視覚情報として得ることができる。近年、マルチスライスＣＴによる高速スキャンが可能になり益々造影ＣＴのダイナミックスキャンの活用範囲が拡大していくものと考えられる。
【０００３】
その一つの方向性として、脳組織内の毛細血管の血流動態に関するインデックスを演算するためのＣＢＰスタディと呼ばれる方法がある。ＣＢＰスタディでは、組織内の局所的な血流動態、つまり局所組織内の毛細血管を通過する血流の動態を定量的に表すＣＢＰ、ＣＢＶ、ＭＴＴ、Ｅｒｒ等のインデックスを求め、またこれらインデックスのマップを出力する。
【０００４】
ＣＢＰは、脳組織の毛細血管内の単位体積及び単位時間あたりの血流量［ml/１００ml/miｎ］を表し、ＣＢＶは、脳組織内の単位体積あたりの血液量［ml/１００ml］、ＭＴＴは毛細血管の血液平均通過時間［秒］を表し、Ｅｒｒは伝達関数を近似する際の残差の総和又は残差の２乗和の平方根を表している。
【０００５】
これら脳組織中の毛細血管の血流動態を定量的に表しているインデックスＣＢＰ、ＣＢＶ、ＭＴＴは、脳虚血卒中が発症してからの経過時間情報ととともに、虚血性脳血管障害の病体鑑別、毛細血管の拡大の有無、血流速などの評価のための有益な情報として期待されている。例えば、一般に虚血性の脳血管障害では、提供する動脈の血圧が低下し、その血管内の血流速の低下が見られる。その結果、ＣＢＶは一定でも、ＭＴＴが延長し、ＣＢＰは低下する。また、脳梗塞超急性期では、血圧低下による血流速の低下を補うために、毛細血管を拡張させ、血流速を増加させることにより、血流量ＣＢＰの低下を抑制しようとする働き（オートレギュレーション）がある。従って、ＭＴＴが延長することにより、ＣＢＰが低下しても、ＣＢＶが増加していれば、毛細血管の再開通の可能性を示唆する情報となる。
【０００６】
ＣＢＰスタディではトレーサーとして脳血管透過性を持たない造影剤、例えばヨード造影剤が使用される。ヨード造影剤は例えばインジェクターにより肘静脈から注入される。インジェクターにより静注されたヨード造影剤は、心臓、肺を経由して、脳動脈へ流れ込む。そして、造影剤は、脳動脈から、脳組織内の毛細血管を経て、脳静脈へと流れ出ていく。このとき、ヨード造影剤は正常な脳組織内の毛細血管では血管外へ漏れ出ることなく通過する。図１はこの様子を模式的に示している。
【０００７】
造影剤の通過の様子をダイナミックＣＴで撮影して、その連続画像から、脳動脈上の画素の時間濃度曲線Ｃａ（ｔ）、脳組織（毛細血管）上の画素の時間濃度曲線Ｃｉ（ｔ）、脳静脈上の画素の時間濃度曲線Ｃｓｓｓ（ｔ）をそれぞれ測定する。
【０００８】
ここで、ＣＢＰスタディでは、脳動脈の濃度時間曲線Ｃａ（ｔ）と脳組織の濃度時間曲線Ｃｉ（ｔ）との間で成り立つ理想的な関係を解析モデルとしている。脳組織に入る直前の血管から造影剤を注入した場合、脳組織の単位体積（１画素）内の時間濃度曲線は立ち上がりが垂直で、しばらくは一定の値を維持し、その後、急勾配で立ち下がる形になる。これを矩形関数で近似する（ｂｏｘ−ＭＴＦ法：ｂｏｘ−ＭｏｄｕｌａｔｉｏｎＴｒａｎｓｆｅｒＦｕｎｃｔｉｏｎｍｅｔｈｏｄ）。
【０００９】
つまり、脳動脈の時間濃度曲線Ｃａ（ｔ）を入力関数、脳組織の時間濃度曲線Ｃｉ（ｔ）を出力関数として、入力関数と出力関数との間の伝達関数を矩形関数で近似する。伝達関数は、トレーサーが毛細血管を通過する過程を表している。
【００１０】
ＣＢＰスタディの問題は次の通りである。
上記ＣＢＰ、ＣＢＶ、ＭＴＴ、Ｅｒｒの各インデックスは、各画素（x,y,z）ごとに算出されるので、その値を画素値とする画像を構成することができ、このような画像をマップと呼ぶ。たとえばＲ種のインデックスが得られる場合、Ｒ枚のマップが構成できる。このようにして作成されたＲ枚のマップは、各画素がベクトル値をとる１枚のマップ（ベクトル値マップ）と見なすことができる。すなわち、次のように表せる。
Ｖ_ｋ（x,y,z）＝<Ｐ_ｋ _, _１（x,y,z）, Ｐ_ｋ _, _２（x,y,z）, ... , Ｐ_ｋ _, _Ｒ（x,y,z）>
例えばＣＢＰスタディでは、典型的には、Ｒ＝４とし、Ｐ_ｋ _, _１（x,y,z）はＣＢＰの値を、Ｐ_ｋ _, _２（x,y,z）はＣＢＶの値を、Ｐ_ｋ _, _３（x,y,z）はＭＴＴの値を、Ｐ_ｋ _, _４（x,y,z）は残差Ｅｒｒの値を表すように構成できる。
【００１１】
このようなベクトル値マップＶ_ｋは、参照した脳動脈の時間濃度曲線Ｃａ（ｔ）_ｋごとに作成される。例えば、左右の中大脳動脈、前大脳動脈、後大脳動脈から脳動脈の時間濃度曲線を得たとするとＫ＝６、さらに病変部周辺にある動脈数カ所から脳動脈の時間濃度曲線を得たとすると、Ｋ＝１０〜１５程度になる。
【００１２】
このように脳動脈の時間濃度曲線Ｃａ（ｔ）_ｋ（ｋ＝１,２，・・・，Ｋ）の数Ｋが大きい場合に、結果として得られるベクトル値マップＶ_ｋ（ｋ＝１,２，・・・，Ｋ）の枚数が多いために、観察するのに不便である。すなわち通常のグレースケール画像あるいはカラースケール画像として観察しようとすれば、一つのマップがＲ枚の画像から構成され、これがＫ組あるのだから、合計Ｋ×Ｒ枚の画像を比較しなくてはならない。さらに、どの部位がどの動脈によって栄養されているのかは必ずしも自明でなく、解剖学的知識を用いて、各部位ごとにどのマップＶ_ｋ（ｋ＝１,２，・・・，Ｋ）を観察するべきかを判断しなくてはならない。特に脳梗塞等の脳血管障害を生じている症例においては、組織を支配（depend）しているのがどの動脈かは、必ずしも解剖学的知識とは一致せず、異常な支配がしばしば見られる。これらの問題によって、ベクトル値マップの読影が難しいという問題点がある。
【００１３】
また、マルチスライスＣＴ、あるいはボリュームＣＴによって撮影されたダイナミックＣＴ画像においては、さらに多数の動脈が観察される。これは同じ動脈が複数のスライス中で観察できるからである。これらの動脈の断層像全部について脳動脈の時間濃度曲線を作ると非常に数が多くなる。
【００１４】
また、ＣＢＰスタディには次のような問題もある。肘静脈に、ボーラスインジェクションを行った場合、ＣＴで観察される造影効果は、血液のＣＴ値（造影されない時数十ＨＵ）が最大数百ＨＵに上昇する。しかし、脳血流を有効に解析するためには造影効果を高々数パーセント以内の誤差で計測できなくてはならない。すなわち、血液の造影効果（ＣＴ値の上昇）が２０〜４０ＨＵ程度であってもこれを検出できる必要がある。
【００１５】
単位体積の脳組織中に占める毛細血管の体積比率は高々３〜４パーセント程度である。従って、血液のＣＴ値が２０〜４０ＨＵ上昇した場合、脳組織の平均ＣＴ値は、０．５〜１．５ＨＵ程度上昇するに過ぎない。
【００１６】
ＣＴ画像ではノイズの標準偏差（ｓｄ）は、照射Ｘ線量の平方根に反比例し、典型的な照射条件においてｓｄは例えば５〜１０ＨＵ程度である。従って、０．５ＨＵの造影効果を検出するためには、Ｘ線量を１０〜１００倍程度増やさねばならず、これは患者の被曝線量が著しく大きくなることを意味する。また、ダイナミックＣＴにおいては同一箇所を数十回に渡って撮影するのであるから、撮影箇所に於ける皮膚の被曝は通常の数百〜数千倍に至ることになり、炎症・脱毛・壊死・発癌等の放射線障害を考慮すると、現実的ではない。
【００１７】
むしろダイナミックＣＴにおいてはＸ線量を通常の撮影よりも減らさなくてはならない。一般に、１スキャン当たりのＸ線量を例えば通常の１／２〜１／１０程度に減じることが行われる。これによって、通常の１回のＣＴ撮影に比べて数倍〜２０倍程度のＸ線被曝に留めることができ、これは放射線障害を生じない程度である。しかし、このようなＸ線量を低減したＣＴ画像において、ｓｄは例えば１５〜２０ＨＵ程度であり、０．５〜１．５ＨＵ程度の造影効果は到底検出できない。
【００１８】
そこで、画像のノイズ成分を抑制することが、ＣＢＰスタディでは重要な課題の１つである。そのため、１）スライス厚を厚くする、２）近隣画素を平均化する、３）画像を平滑化処理に通す、が一般的に取りうる方策である。しかしこれらには以下のような問題点がある。
【００１９】
“スライス厚を厚くする”ために、撮影時にスライス厚を厚く設定するか、連続する薄いスライスの画像数枚を平均して厚いスライスの画像を生成する。スライス厚に比例して画素当たりのＸ線量が増えるため、画像ノイズのｓｄは、スライス厚の平方根に反比例して小さくなる。しかしながら、スライス厚を厚くすることによって、パーシャルボリューム効果が生じ、すなわち１個の画素が、一様な脳組織を表しておらず、複数の組織（白質、灰白質、血管、脳溝、脳室など）の平均的なＣＴ値を表すことになる確率が大きくなり、解析結果として得られる脳血流量等の値の誤差が大きくなる。
【００２０】
特に血管の影響を含む画素は、正常な解析が不可能である。このためスライス厚を厚くすると、不正確で、しかも解析不可能な画素を沢山含む非常に品質の悪い結果しか得られなくなる。
【００２１】
“近隣画素を平均化する”では、空間解像度が或る程度犠牲になる。例えば一辺がｎ個の画素からなる正方形の領域（ｎ×ｎ個の画素を含む）の平均値を求め、これをその正方形全体の平均ＣＴ値とし、このような正方形を画素とみなし、これを敷き詰めて「画素束ね画像」を構成する。もとの画像が例えば一辺５１２個の画素からなる（５１２×５１２個の画素を含む）とし、ｎ＝２とすれば、「画素束ね画像」は一辺（５１２／２）個の画素から構成される（２５６×２５６個の画素を含む）画像となる。この方法によれば、ノイズは、ｎに反比例して減少させることが可能である。さらに、解析対象となる画素の数が１／（ｎ×ｎ）倍になるため、計算量も小さくなるという利点がある。
【００２２】
しかしながら、ｎを大きくすると、空間解像度が低下し、それに伴ってパーシャルボリューム効果が生じ、すなわち１個の画素が、一様な脳組織を表して居らず、複数の組織（白質、灰白質、血管、脳溝、脳室など）の平均的なＣＴ値を表すことになる確率が大きくなり、解析結果として得られる脳血流量等の値の誤差が大きくなる。特に血管の影響を含む画素は、正常な解析が不可能である。このため、ｎを大きくすると、空間解像度が低く、不正確で、しかも解析不可能な画素を沢山含む非常に品質の悪い結果しか得られなくなる。このため、実用上は、ｎ＝２〜４程度が限界であり、これだけでは十分なノイズ抑制効果が得られない。
【００２３】
また、画像の平滑化、すなわち１枚のＣＴ画像ごとに、２次元の空間フィルタを作用させて平滑化を行う方法を用いると、十分なノイズ抑制効果と引き換えに、空間解像度が著しく損なわれる。特に、太い血管（動脈・静脈）が存在する箇所に近接している画素には、太い血管において生じた造影効果の影響が及ぶことになり、これらの画素の時間濃度曲線は正しくなくなってしまう。従ってごく軽度の平滑化を行うに留めねばならない。ここで、ごく軽度の平滑化を行うに際して重要なのは、画像フィルタのサイズをごく小さくする事、例えば、３×３程度に設定することである。３×３の平滑化フィルタを用いて最大の画像ノイズ抑制効果を得ようとすると、その上限は、ノイズｓｄを１／３に低減することであり、それ以上にノイズを抑制するのは不可能である。従って十分なノイズ抑制効果は得られない。
【００２４】
一方、時間的平滑化、すなわち各画素について得られた時間濃度曲線を曲線とみなして、これを１次元フィルターで平滑化する手法を用いると、十分なノイズ抑制効果を得ようとすると時間分解能を著しく損なう。元来、ＣＢＰスタディでダイナミックＣＴを行うのは短いサンプリング周期で撮影を行うことによって高い時間分解能を得て、時間濃度曲線の僅かで速い変化（特に生理学的構造に起因する平滑化効果がどの程度生じているか）を精密に計測することが目的であるから、時間的平滑化は全く適当でない。
【００２５】
【発明が解決しようとする課題】
本発明の目的は、ＣＢＰスタディにより得られるマップの診断効率を向上することにある。
本発明の他の目的は、空間及び時間分解能の低下を抑えて、ノイズを抑制することにより、ＣＢＰスタディの解析精度を向上することにある。
【００２６】
【課題を解決するための手段】
本発明のある局面は、造影剤を注入された被検体の特定部位に関する連続的な複数の画像から、前記特定部位内の動脈の設定を支援するために、前記画像を構成する複数の画素各々の時間濃度曲線の特性値を計算し、特性値マップを表示し、前記動脈に関する第１時間濃度曲線と、前記特定部位内の組織に関する第２時間濃度曲線を作成し、前記動脈各々に対する前記組織の局所血流動態を表す伝達関数を、前記伝達関数と第１時間濃度曲線とのコンボルーションに対する前記第２時間濃度曲線の差が最小化するように計算し、前記伝達関数から、前記動脈に対する前記局所血流動態に関するインデックスを計算して表示することを特徴とするインデックス計算方法を提供する。
【００２７】
【発明の実施の形態】
以下、図面を参照して本発明を好ましい実施形態により説明する。
本実施形態の特徴としては、ＣＢＰスタディにより生成される多くのインデックスマップを１枚に合成することにより、ＣＢＰスタディの診断効率を向上すること、さらにコヒーレントフィルタを用いることにより、ノイズの低減と、空間及び時間分解能の低下抑制とを両立することによりインデックスの精度を向上することにある。
【００２８】
なお、本実施形態は、被検体の特定部位に関する時間的に連続する複数枚の画像から、局所の血流動態を表すインデックスを計算する方法及び装置に関するものであり、対象となる複数の画像を発生するモダリティは、特定装置に限定されることは無く、例えば、Ｘ線コンピュータトモグラフィ装置（Ｘ線ＣＴ装置）、シングルフォトンエミッショントモグラフィ装置（ＳＰＥＣＴ）、ポジトロンエミッショントモグラフィ装置（ＰＥＴ）、磁気共鳴イメージング装置（ＭＲＩ）のいずれでも良い。なお、ここでは、Ｘ線ＣＴ装置を例に説明する。
【００２９】
（装置構成）
図２には、本実施形態に係るＸ線ＣＴ装置の構成を示している。Ｘ線ＣＴ装置は、ガントリ部１０とコンピュータ装置２０とから構成される。ガントリ部１０は、Ｘ線管１０１、高電圧発生装置１０１ａ、Ｘ線検出器１０２、データ収集部１０３（ＤＡＳ；ＤａｔａＡｑｕｉｓｉｔｉｏｎＳｙｓｔｅｍ）とを有する。Ｘ線管１０１とＸ線検出器１０２とは、高速で且つ連続的に回転する図示しない回転リングに被検体Ｐを挟んで互いに対向する位置に搭載される。
【００３０】
コンピュータ装置２０は、画像処理装置３０と、画像表示部１０７と、入力部１０９とから構成される。画像処理装置３０は、制御部１０８を中枢として、データ収集部１０３から出力される生データを補正処理等を経て投影データに変換する前処理部１０４、投影データを記憶するメモリ部１０５、投影データからＣＴ画像データを再構成する画像再構成部１０６、ＣＴ画像データを保管する記憶装置１０Ｍ、ＣＴ画像データに対してコヒーレントフィルタ処理を実行するコヒーレントフィルタ処理部１１０、及びコヒーレントフィルタ処理を受けたＣＴ画像データを使ってＣＢＰスタディ処理を実行するＣＢＰスタディ処理部１２０とから構成される。
【００３１】
コヒーレントフィルタ処理部１１０は、分散値推定部１１１、重み関数演算部１１２、画素値演算部（コヒーレントフィルタ部）１１３とから構成される。これら分散値推定部１１１、重み関数演算部１１２、画素値演算部１１３の機能については後述するコヒーレントフィルタ処理の詳細説明の中で説明する。
【００３２】
ＣＢＰスタディ処理部１２０は、ＲＯＩ設定支援部１２１、時間濃度曲線作成部１２２、脳動脈時間濃度曲線補正部１２３、ＭＴＦ処理部１２４、インデックス計算部１２５、マップ作成部１２６、マップ合成部１２７から構成される。
【００３３】
ＲＯＩ設定支援部１２１は、ＣＴ画像上に脳動脈や脳静脈に対して関心領域ＲＯＩを設定する作業を支援するための情報（脳動脈ＲＯＩのためのＡＴマップ、ＰＴマップ、ＴＴマップ等）を作成し提供する。
【００３４】
なお、脳動脈ＲＯＩは、例えば前大脳動脈（ＡＣＡ）、中大脳動脈（ＭＣＡ）、後大脳動脈（ＰＣＡ）を対象として、左脳、右脳それぞれの領域に個別に設定される。従って、この例では左右に３個ずつ、合計６個の脳動脈ＲＯＩが設定される。また、脳動脈の時間濃度曲線Ｃａ（ｔ）を補正するために、他の時間濃度曲線Ｃｓｓｓ（ｔ）が利用される。この時間濃度曲線Ｃｓｓｓ（ｔ）は、パーシャルボリュームを含まない画素が存在するのに充分に太い血管上に設定されたＲＯＩに関して作成される。Ｃｓｓｓ（ｔ）のＲＯＩは、例えば、脳血管の中で最も太い上矢状静脈洞に設定される。
【００３５】
時間濃度曲線作成部１２２は、記憶装置１０Ｍに記憶されているダイナミックＣＴ画像データ（時間的に連続した複数枚の画像データ）から脳動脈、脳静脈及び脳組織（毛細血管）に関する時間濃度曲線を作成する。なお、脳動脈の時間濃度曲線Ｃａ（ｔ）は、設定された例えば６つの脳動脈ＲＯＩに関して個々に作成される。脳静脈の時間濃度曲線Ｃｓｓｓ（ｔ）は、上矢状静脈洞に設定された脳静脈ＲＯＩに関して作成される。また、脳組織の時間濃度曲線Ｃｉ（ｔ）は、脳組織上の全画素を対象として画素ごとに作成される。
【００３６】
脳動脈時間濃度曲線補正部１２３は、ノイズやパーシャルボリューム効果の影響を除去するために、脳動脈の時間濃度曲線Ｃａ（ｔ）を、上矢状静脈洞の時間濃度曲線Ｃｓｓｓ（ｔ）に基づいて補正する。この補正方法については後述する。ＭＴＦ処理部１２４は、補正された脳動脈の時間濃度曲線Ｃａ（ｔ）と、脳組織の時間濃度曲線Ｃｉ（ｔ）とに基づいて、ｂｏｘ−ＭＴＦ法により、伝達関数ＭＴＦを、脳組織領域内の全画素を対象として画素ごとに計算する。
【００３７】
インデックス計算部１２５は、計算された伝達関数ＭＴＦから脳組織の血流動態を表すインデックス（ＣＢＰ、ＣＢＶ、ＭＴＴ、Ｅｒｒ）を、脳組織領域内の全画素を対象として画素ごと計算する。マップ作成部１２６は、計算されたインデックス各々のマップを、脳動脈（ＡＣＡ，ＭＣＡ，ＰＣＡ）ごとに生成する。マップは、各スライスに関して、インデックスの種類（＝４）×脳動脈の数（ＡＣＡ，ＭＣＡ，ＰＣＡの３つ）＝１２種類作成される。マルチスライスでは、そのスライス数倍の種類のマップが作成される。この膨大な枚数のマップの枚数を合成処理により減らして診断効率を向上させるためにマップ合成部１２７が設けられている。
【００３８】
以下に、コヒーレントフィルタ処理とＣＢＰスタディ処理について順番に説明する。
コヒーレントフィルタの原理について簡単に説明すると、近傍の例えば３×３等の局所内画素を加重平均し、その加重平均値を局所中心画素の値とすることを基本として、周辺画素各々の重みを中心画素と周辺画素との間の類似度に従って変えることを特徴としたものである。ここで言う類似度とは、画素間で、解剖学的に近い組織、具体的には同じ脳動脈の支配下にある脳組織（毛細血管）どうしである可能性の度合いを示す指標であり、この類似度が高い画素に対しては高い重みを与え、逆に類似度が低い画素に対してはゼロに近い低い重みを与えることにより、ノイズ抑制を果たしながらも、空間分解能の低下を抑制することを可能としている。なお、以下では、類似度は、適宜、適合度又は危険率に言い換えられる。
【００３９】
本実施形態では、脳血管透過性を持たない造影剤、例えばヨード造影剤を注入（静注）された被検体の脳を撮影対象として連続的に取得した複数のＣＴ画像（ダイナミックＣＴ画像）を用いて、各画素の時間濃度曲線の比較により類似度を計算する。そのため類似度の確からしさは、サンプリング周波数、つまり単位時間あたりの画像枚数と、サンプリング数、つまり全画像枚数とに依存して決まる。そこでスキャン間隔を例えば０．５秒に短縮することが効果的である。
【００４０】
（コヒーレントフィルタ）
（コヒーレントフィルタの一般的説明）
（画素値ｖ（ｘ））
一般に、カメラやＣＴスキャナ等の撮像手段を介して取得されたデジタル画像は、複数の画素（ｐｉｘｅｌ）から構成されている（あるいは、当該画像をそのような画素の集合として考えることができる。）。以下の説明では、当該画素の位置をベクトルｘ（すなわち座標値のベクトル）として表し、画素ｘが有する値（例えば濃淡を表わす数値、ＣＴ値ＨＵ）をＫ次元ベクトルとして表す。２次元画像の場合、画素ｘとは画像上における位置を表す座標値（ｘ、ｙ）を示す２次元ベクトルである。ある画素ｘについて定義される「画素値ｖ（ｘ）」を、
ｖ（ｘ）＝（ｖ_１（ｘ），ｖ_２（ｘ），…，ｖ_Ｋ（ｘ）） … （１）
と表記する。なお、この（１）式の右辺における、ｖ_１（ｘ），ｖ_２（ｘ），…，ｖ_Ｋ（ｘ）それぞれを、以下では、画素ｘについての「スカラー値」と呼ぶことにする。
【００４１】
例えば、画像が「カラー画像」であるとき、各画素が、それぞれ三原色（赤，緑，青）の明るさ（スカラー値）を有することから、これら各画素の画素値ｖ（ｘ）は、その次元がＫ＝３のベクトルであると考えることができる。すなわち上記（１）式の右辺各項添え字が例えば「赤」，「緑」及び「青」に対応付けられる。また例えば、画像がＫ枚の静止画像から構成される動画像であって、第ｎ番目の画像の各画素はスカラー値ｖ_ｎ（ｘ）を持つという場合には、Ｋ枚の静止画像上、共通の同一点（同一座標）の画素ｘの持つ画素値（スカラー値）を並べて構成される、Ｋ次元ベクトル値ｖ_ｎ（ｘ）＝（ｖ_１（ｘ），ｖ_２（ｘ），…，ｖ_Ｋ（ｘ）が以下で述べるベクトル値としての画素値である。
【００４２】
（類似度（適合度又は危険率）と重み）
上記画素ｘに対して、適当な周辺画素の集合Ｎ（ｘ）を考える（この集合Ｎ（ｘ）は画素ｘを含んでよい。）。次に、中心画素ｘに対するＮ（ｘ）の要素である周辺画素ｙの重みｗ（ｐ（ｘ，ｙ））を考える。この重みｗ（ｐ（ｘ，ｙ））は、次に記す性質を有する。
【００４３】
（類似度ｐ（ｘ，ｙ））
まず、重みｗ（ｐ（ｘ，ｙ））の値を左右する関数ｐ（ｘ，ｙ）の意味について述べる。このｐ（ｘ，ｙ）は、本実施形態にいう「類似度」を定量化する手段であり、一般的にいえば、中心画素ｘと周辺画素ｙ∈Ｎ（ｘ）とが、何らかの意味でどの程度類似しているか（例えば、両画素ｘ及びｙの上記画素値ｖ（ｘ）及びｖ（ｙ）間に認められる統計的差異の程度）を示す具体的数値を与える。
【００４４】
より具体的には例えば、ｐ（ｘ，ｙ）が大きな値を示すときには、画素ｘと画素ｙとの間に、「統計的に有意な差がなく（つまり類似度が高い）」、類似である可能性が高いと判断され、ｐ（ｘ，ｙ）が小さい値を示すときには、画素ｘと画素ｙとの間に、「統計的に有意な差があり（つまり類似度が低い）」、の如く判断されるということである。
【００４５】
ところで、画素値ｖ（ｘ）及びｖ（ｙ）（ないしスカラー値ｖ_１（ｘ），…，ｖ_Ｋ（ｘ）及びｖ_１（ｙ），…，ｖ_Ｋ（ｙ））には、必ずノイズが含まれている。例えば、画像がＣＣＤ撮像素子により取得された場合を考えると、それを構成する各画素については、素子内の暗電流や外界から入射する光量の不規則変動に起因するノイズ等が存在する。
【００４６】
このようなノイズは、一般に、全画素についてまちまちな値をとるため、画素ｘと画素ｙとが、仮に（外界における）同一物体を反映したものである場合であっても、実際に観測される画像上では、同一の値を持たないことがある。このことを逆にいえば、いずれも同一物体を反映した画素ｘと画素ｙにおいて、それぞれのノイズを除去した状況を仮に想定すれば、これらは該同一物体を表象するものとして画像上に表示され（＝そのように認識され）るし、また、両者は本来同一の（あるいはごく近い）画素値を有する。
【００４７】
そこで、上述したノイズの性質を踏まえ、上記のｐ（ｘ，ｙ）に関し、統計的検定法でよく知られている「帰無仮説」の概念を用いると、このｐ（ｘ，ｙ）については、具体的に次のように言うことができる。すなわち、帰無仮説Ｈ「画素ｘと画素ｙとはそれぞれのノイズを除去した場合に同一の画素値を有する」言いかえれば「ｖ（ｘ）＝ｖ（ｙ）、ただし、両画素のノイズに起因する差異を除く」を立てる（つまり、このような命題が成立する場合、「両画素ｘ及びｙとの間の類似度が高い（適合度が大きい）」と考える。）と、関数ｐ（ｘ，ｙ）は、この仮説Ｈを棄却する場合の危険率（あるいは、有意水準）として構成できる（この場合、ｐ（ｘ，ｙ）は、その値域が[０，１]であるような関数として定義される（ｐ（ｘ，ｙ）∈[０，１]））。
【００４８】
したがって、危険率ｐ（ｘ，ｙ）が大きい場合、すなわち棄却が誤りである危険性が大きい場合には上記仮説Ｈを満たす可能性が高いといえ、逆に小さい場合、すなわち棄却が誤りである危険性が小さい場合には仮説Ｈを満たさない可能性が高いということができる（なお、統計的検定における周知事項ではあるが、仮説Ｈが「棄却」されないといっても、それが「真」であることを意味するわけではない。この場合、仮説Ｈが示す命題が、否定し得ないことを意味するに過ぎない）。
【００４９】
（重みｗ（ｐ（ｘ，ｙ）））
さて、重みｗ（ｐ（ｘ，ｙ））は、その表され方から明らかな通り、上記したような危険率ｐ（ｘ，ｙ）の関数（より一般には、適合度の関数（適合度をρ（ｘ，ｙ）とすれば、ｗ（ρ（ｘ，ｙ））となるように構成できる）であり、また、この重みｗ（ｐ（ｘ，ｙ））を求めるため、ｘ及びｙの組み合わせそれぞれについて求められた危険率ｐ（ｘ，ｙ）に作用させる重み関数ｗは、一般的にいうと、上記「棄却」を具現化する作用を有するものである。具体的には、危険率ｐ（ｘ，ｙ）が大きい場合には重み関数ｗの値、すなわち重みｗ（ｐ（ｘ，ｙ））が大きな正の値をとり、その逆の場合には小さな正の値（又は“０”）をとる、というように調整されている（重み関数ｗの具体的形式については後述する。）。つまり、重みｗ（ｐ（ｘ，ｙ））は、画素ｘと画素ｙとが、上記仮説Ｈに示される命題を満たすらしい場合には、大きい値をとり、その逆の場合には小さい値をとる。一例として特に、ｗのとりうる値が”０”かまたは”０”でない一定値の２通りしかないように構成してもよい。
【００５０】
なお、以上までに述べた仮説Ｈ、危険率ｐ（ｘ，ｙ）、重みｗ（ｐ（ｘ，ｙ））間の関係をまとめると、帰無仮説Ｈが正しい可能性が高いとき、類似度ｐも高くなり、その画素に与える重みｗを高くし、一方、帰無仮説Ｈが正しい可能性が低いとき、類似度ｐも低くなり、その画素に与える重みｗを低くする。このように加重平均値への寄与度（重み）を類似度に応じて変えることにより、分解能の低下を抑えながら、ノイズを効果的に抑制することが可能となる。また、重み関数ｗ（ｔ）は、より一般に、「ｔ∈[０，１］で定義される非負の単調増加関数」ということができ、また、該ｗ（ｔ）の満たすべき性質は、少なくともそのようであればよい。
【００５１】
（コヒーレントフィルタ処理）
以上までの説明により、「コヒーレントフィルタ」は次のように導かれる。すなわちまず、画像を構成するある画素ｘに対し、集合Ｎ（ｘ）の要素たる画素ｙのすべてについて上記した重みｗ（ｐ（ｘ，ｙ））を計算する。次に、これら複数の重みｗ（ｐ（ｘ，ｙ））を用いて、当該画素ｘを構成する新たなスカラー値ｖ´_ｋ（ｘ）を、以下の（２）式で計算する。すなわち、
【数１】

ただし、ｋ＝１，２，…，Ｋである。そして、この式で求められたｖ´_ｋ（ｘ）を用いて、当該画素ｘの変換後の画素値（新たな画素値）ｖ´（ｘ）を、
ｖ´（ｘ）＝（ｖ´_１（ｘ），ｖ´_２（ｘ），…，ｖ´_Ｋ（ｘ）） … （３）
として構成する。
【００５２】
ここに、上記（２）式で表される、画素値ｖ（ｙ）＝（ｖ_１（ｙ），ｖ_２（ｙ），…，ｖ_Ｋ（ｙ））（ｙ＝ｘである場合を含む。）を、ｖ´（ｘ）＝（ｖ´_１（ｘ），ｖ´_２（ｘ），…，ｖ´_Ｋ（ｘ））に変換するフィルタが、「コヒーレントフィルタ」の形式である。これはその表式から明らかな通り、画素値を構成するスカラー値ｖ_ｋ（ｙ）の重み付け平均値を表している。
【００５３】
このような処理は、以下のような結果をもたらす。すなわち、ノイズを除かれた画素ｘの画素値ｖ´（ｘ）と同一の画素値をとることが確からしい（＝上記仮説Ｈの命題を満たす可能性が高い）画素ｙの寄与度を高めた重み付け平均値ｖ´_ｋ（ｘ）から構成されたベクトルを表すこととなる。また、このような画素ｙが十分な数存在するならば、画素値ｖ´（ｘ）は、画素ｘが本来有すべきその真値から外れることなく、上記したような平均化の作用によりノイズを抑制しながらも、分解能は低下しない値を有することとなる。
【００５４】
なお、危険率ｐ（ｘ，ｙ）が小さく、したがって、帰無仮説Ｈが「棄却」され、重みｗ（ｐ（ｘ，ｙ））が小さくなるような場合であっても、上記記述からもわかる通り、必ずしもこれを完全に「棄却」するとは限らない。このようなことは、後述する重み関数ｗの具体的形式に依存するところであるが、危険率ｐ（ｘ，ｙ）が“０”（＝０％）に近いような場合でも、ｗ（ｐ（ｘ，ｙ））≠０（ただし、ｐ（ｘ，ｙ）が“１”に近い場合に比べて、より小さな正の値ではある。）としてよい（なお、ｐ（ｘ，ｙ）＝１である場合とは、後述するように、ｖ（ｘ）＝ｖ（ｙ）のときである。）。
【００５５】
すなわち、完全な棄却ということではなく、小さな寄与は認めるように構成してもよいということである（なおこのような場合に、ｗ（ｐ（ｘ，ｙ））＝０とするのであれば、完全な棄却を行うのと同義である。
【００５６】
このような処理は、一般的に次のように言える。すなわち、ある画像を構成する複数の画素ｘが存在するとき、この画素ｘとある任意の画素ｙ（上記ではｙ∈Ｎ（ｘ）とされた）との適合度ｐ（ｘ，ｙ）を定量化し、該適合度が大きい場合には、画素値ｖ（ｙ）を利用した重み付き平均化処理において、当該画素ｙの大きな寄与を認め、適合度が小さい場合には小さな寄与しか認めないようにすることで、当該画素ｘのノイズを有効に抑制する画像処理方法である、といえる。いわば、画素ｘと画素ｙとが「似たもの同士」のときには、該画素ｙを前記平均化処理に、より貢献させ、「似ていないもの同士」のときには、該画素ｙを全く又は殆ど無視する（重みをゼロ又はその近似値）、と言い換えてもよい。
【００５７】
このような処理を画像全体に施すことにより、画像のぼけ、つまり空間的分解能の低下を殆ど生じることなく、極めて高いノイズ抑制効果を発揮することができる。また、ノイズ抑制という用途に限定せず、例えばパターン認識の分野においても、重み関数、あるいはコヒーレントフィルタを好適な具体的形式にすることによって、優れた効果を発揮することができる。
【００５８】
ここで上記した「ダイナミックＣＴ」撮影とは、上記Ｘ線管１０１及びＸ線検出器１０２が被検体Ｐの同一部位を反復撮影（反復スキャン、連続回転型ＣＴ装置では、連続回転による反復撮影がしばしば行われる。）して次々に投影データを取得するとともに、該投影データに基づいて次々に再構成処理を行って時系列的な一連の画像を得る撮影方式のことをいう（この場合、画像表示部１０７における画像表示は、例えば図示しないカウンタ等によって、その画像の元となった投影データ収集に係るスキャン開始点又は終点から一定時間後に行われるように制御される）。
【００５９】
したがって、このように取得・表示される画像は、映画等と同様に時系列的な複数枚の静止画像からなる、いわゆる動画像となる。なお、このような撮影方式は、典型的には、被検体Ｐに対し造影剤を注入し、その経時変化を観察・解析して、例えば血管における狭窄や閉塞等その他病変部の病態を分析するために用いられる。また、造影剤投与の前後２回だけに限り同一部位のＣＴ撮影を行う方式も、広義のダイナミックＣＴ撮影と考えることができる。
【００６０】
さて、従来においては、上記のような「ダイナミックＣＴ」撮影時、例えばＫ回の撮影を実施する間に被検体Ｐに何らかの変化（例えば、造影剤の濃度変化や呼吸動等が一般的に考えられる）があった場合、空間解像度を損なわず画像ノイズを抑制するためには、時間方向の平滑化を行うほかなかった。その結果、時間分解能が損なわれるという弊害は避け得なかった。
【００６１】
ところが、ダイナミックＣＴ撮影により取得される画像は、上述したように動画像であって、時間的変化を仔細に観察する目的で撮影を行うものであるから、その時間分解能が損なわれるというのは、本来、好ましい状況とは言えない。
【００６２】
コヒーレントフィルタを利用すれば、時間分解能を損ねず、Ｋ枚の静止画像のすべて（複数枚の画像）からノイズを抑制することが可能な、次のようなダイナミック・コヒーレントフィルタ処理を実施することができる。
【００６３】
まず、上記のようにして得られた動画像たるＫ枚の静止画像につき定義される画素ｘについては、既に述べたように、画素値ｖ（ｘ）として、
ｖ（ｘ）＝（ｖ_１（ｘ），ｖ_２（ｘ），…，ｖ_Ｋ（ｘ）） … （１再掲）
を構成することができる。ここで右辺各項における添え字１，２，…，Ｋは、Ｋ枚の各静止画像の通し番号である。
【００６４】
次に、この場合における重み関数ｗ１の具体的形式を、例えば次の（４）式により与える。
【００６５】
【数２】

【００６６】
ただし、ｙ∈Ｎ（ｘ）であって、かつ、この集合Ｎ（ｘ）は、画素ｘにつき任意に設定してよい（＝どのような基準によって設定してもよい。）。しかし実際上は、画素ｘと該画素ｘから遠く離れた位置にある画素ｙとが仮説「ｖ（ｘ）＝ｖ（ｙ）。ただし、両画素のノイズに起因する差異を除く」を満たす可能性は一般に低いといえるから、集合Ｎ（ｘ）をｘに近接している画素の集合という基準で限定することは、演算速度向上等の実用的な意義がある。
【００６７】
したがってここでは、その一例として、集合Ｎ（ｘ）を、当該画素ｘを中心としたその周囲の矩形状エリアに含まれる画素の集合、とする。より具体的に、集合Ｎ（ｘ）としては、例えば、いま注目している静止画像一枚を構成する全画素が１２８×１２８画素であるような場合に、前記画素ｘを中心とした３×３画素分のエリアとしたり、また、５１２×５１２画素であるような場合に、当該画素ｘを中心とした１３×１３画素分のエリア等としてもよい。
【００６８】
また、上記（４）式におけるσ_ｋは、ｋ枚目の静止画像の各画素が、そのどれにも共通な一定の程度で有するものと仮定して推定されたノイズの標準偏差であり、一方Ｃは、重みｗ１（ｐ（ｘ，ｙ））が、上記（４）式に代入された場合における作用の程度を決定する調節可能なパラメータである。
【００６９】
以下、これらσ_ｋ及びＣについての説明を順に行う。
まず、（４）式におけるσ_ｋについて説明する（以下では、分散σ_ｋ ^２として説明する）。このσ_ｋ ^２は、上述したように、ｋ枚目の静止画像上の各画素のスカラー値が有するノイズ成分の分散である。そしてまた、上記（４）式における分散σ_ｋ ^２は、ｋ枚目の画像の各画素のスカラー値について一定値たる分散σ_ｋ ^２を持つノイズを含んでいるものと仮定して推定したものである。一般に、このような仮定は、次に記すようなことを背景として、十分な正当性を持つ。
【００７０】
被検体Ｐの大きさ、Ｘ線管１０１及びＸ線検出器１０２、再構成部１０６等の構造が一定で、かつ、照射Ｘ線のエネルギを一定にした状態では、ＣＴ画像のノイズは、照射Ｘ線量、すなわちこれと比例関係にあるＸ線管１０１における管電流と照射時間との積（いわゆる管電流時間積（ｍＡ・ｓ））によって決定される。
【００７１】
一方、ＣＴ画像のノイズは加法的であり、概ねガウス分布に従うことも知られている。すなわち、ある画素ｘの画素値ｖ（ｘ）を構成する任意のスカラー値ｖ_ｎ（ｘ）（ｎ＝１，２，…，Ｋ）について、その真値（ノイズの寄与分を除去した値）をｖ_ｎ ^０（ｘ）とすると、これらの差の値ｖ_ｎ（ｘ）−ｖ_ｎ ^０（ｘ）は、概ね平均０、分散σ_ｋ ^２のガウス分布に従う（なお、照射Ｘ線量ないし管電流時間積ｍＡ・ｓとノイズの分散σ_ｋ ^２とは、概ね反比例関係にある。）。
【００７２】
また、この分散σ_ｋ ^２は、画素ｘの位置そのもの（上で述べたように、例えば各座標値ｘ＝（ｘ，ｙ））にも依存するが、通常のＸ線ＣＴ装置１００においては、Ｘ線管１０１及びＸ線検出器１０２の間に、Ｘ線照射量を調節する物理的なＸ線フィルタ（例えば銅箔や金属塊等により構成された、いわゆる「ウェッジ」あるいは「Ｘ線フィルタ」と呼称されるもの）を備えているため、これを無視することができる。なぜならばウェッジは、被検体Ｐが水とほぼ同じ密度を持つ物質から構成されていることを利用して、どのＸ線検出器１０２においても同程度のＸ線量が検出されるよう、照射されるＸ線の一部を減弱する作用を有するものである。従ってこのようなウェッジによれば、結果的に、ノイズの分散σ_ｋ ^２を画素ｘの位置に殆ど依らない概ね一定値にする効果を生じるからである（ちなみに、このウェッジは、一般に、Ｘ線検出器１０２のダイナミックレンジを有効に利用することを本来の目的として設置されるものである）。
【００７３】
以上のことから、ダイナミックＣＴ撮影により取得されたＫ枚の静止画像上においては、ｋ枚目の静止画像上におけるすべての画素について、分散σ_ｋ ^２がほぼ一定であると推定することは妥当である。むろん、画素ごとに分散が異なる場合について本実施形態を拡張することも容易に推考できる。
【００７４】
さて次に、上記（４）式を具体的に演算するためには、その分散σ_ｋ ^２として、どのような数値をあてるか、が問題となる。このようなことが問題となるのは、通常、ノイズの分布の形は想定できても（上記ではガウス分布）、分散σ_ｋ ^２の具体値は不明であることが多いからである。
【００７５】
更に、一般的に、毎回の撮影毎に照射線量（Ｘ線管電流×照射時間（ｍＡ・ｓ））を変更して撮影を行ってもよく、この場合、画像毎に分散σ_ｋ ^２は異なる。
【００７６】
さて、ｋ枚目の画像（ｋ＝１，２，…，Ｋ）に於いて各画素のスカラー値が持つノイズの分散をσ_ｋ ^２とし、ｋ枚目の画像の撮影に用いた照射線量をＲ_ｋとするとき、σ_ｋ ^２はＲ_ｋに比例する。従って少なくともひとつのｋ＝ｋ_０についてσｋ_O ^２が指定できれば、他のｋに関しても、
【数３】

によってσ_ｋ ^２を正確に推定することができる。
【００７７】
本実施形態に於いては少なくともひとつのｋについて、以下のような方法でσ_ｋ ^２の具体的数値の推定を行うことができる。
【００７８】
Ｋ回の撮影のうち、被検体Ｐに殆ど変化がなかったと仮定することのできるＮ回（１＜Ｎ≦Ｋ）の画像を用いて、実測により、分散σ_ｋ ^２に対する期待値Ｅ[σ_ｋ ^２]を求める方法が有効である。以下説明を簡単にするために、これらＮ枚の画像における照射線量は同じであり、従ってｋ＝１，２，…Ｎに関してσ_ｋ ^２は一定（σ^２と書く）と仮定する。これらＮ枚の画像における、ある画素ｘ_ｆの画素値ｖ（ｘ_ｆ）を構成する各スカラー値ｖ_１（ｘ_ｆ），ｖ_２（ｘ_ｆ），…，ｖ_Ｋ（ｘ_ｆ）が含むノイズは、上述したように平均０、分散σ^２のガウス分布に従うと予想されるから、これらの平均値を以下の（６）式、
【数４】

を用いると、真の分散σ^２に対する期待値Ｅ[σ^２]を、
【数５】

として求めることができる。そして、この分散の期待値Ｅ[σ^２]は、上述した通り、Ｋ枚すべての静止画像上の全画素ｘにつき妥当するものと考えることができ、真の分散σ^２の代用として用いるのに、一定程度以上確からしさが保証された値である。したがって、上記（４）式の実際の演算においては、このＥ[σ^２]を（４）式のσ^２に代入すればよい。
【００７９】
なお、このようなＥ[σ^２]は、より具体的には、Ｋ枚の静止画像中、例えば１枚目と２枚目の静止画像に基づく実測値により求めてもよい（上記（６）及び（７）式で言えば、Ｎ＝２とすることに該当する。）。また、上記（６）及び（７）式の実際の演算に供される画素ｘ_ｆについては、例えば、空気や骨が撮像されている部分を除いた適当な画素ｘ_ｆのみを選定する（複数選定した場合は得られるＥ[σ^２]すべての平均をとる）等といった工夫を施してもよい。さらに、その他一般的には、被検体Ｐの動きによる影響を抑える工夫等を施すと尚よい。
【００８０】
これらＮ枚の画像の撮影において照射線量が一定でない場合においても、σ_ｋ ^２がＲ_ｋに比例することを利用して正しくσ_ｋ ^２を推定することは容易に推考できるであろう。
【００８１】
さて次に、上記（４）式におけるパラメータＣについての説明を行う。まず、（４）式においては、上記一般的形態で述べた危険率ｐ（ｘ，ｙ）の考え方が、以下のように含まれている。すなわち、（４）式の右辺分子における根号内の表式は、いわゆるχ二乗分布に従うとされる当該χ^２値に一致するものであり、これを（２σ）^２で除し、括弧の全体をｅの肩に置いた値は、危険率ｐ１（ｘ，ｙ）そのものである。つまり、
【数６】

そして、上記（４）式は、この（８）式のように表されるｐ１（ｘ，ｙ）に関し、
【数７】

としたものに他ならない。尚、Ａは定数でｐ１が（０〜１）の値になるように規格化されたものである。
【００８２】
結局、（４）式においては、上記したような一般的形態で述べた危険率ｐ（ｘ，ｙ）が陽には表示されてはいないが、重みｗ１（ｐ（ｘ，ｙ））の実態は、上述したように、まさしく危険率（＝ｐ１（ｘ，ｙ））の関数であると見ることができ（（９）式）、すなわち「適合度の関数」である（ただし、危険率と適合度とは、上述したように、一方が増えれば他方も増加する関係にある）。
【００８３】
そして、上記（９）式からわかるように、パラメータＣは、重みｗ１（ｐ（ｘ，ｙ））が、危険率ｐ１（ｘ，ｙ）にどの程度敏感に反応するかを決める効果がある。つまり、Ｃを大きくすると、ｐ１（ｘ，ｙ）がわずかに小さくなるだけで、ｗ１（ｐ（ｘ，ｙ））は０に近づく。また、Ｃを小さくするとそのような過敏な反応を抑制することができる。なお、Ｃとして、具体的には１乃至１０程度とすればよく、好適にはＣ＝３とするとよい。
【００８４】
本実施形態においては、中心画素ｘと周辺画素ｙとの間の類似判定、言い換えると、両画素ｘ及びｙに関する上述した帰無仮説Ｈの棄却の判定は、上述したことから明らかなように、上記危険率ｐ１（ｘ，ｙ）に基づいて、いわゆるχ二乗検定法（統計的検定法）によって決定されている。
【００８５】
また、上記（４）式の表式からわかるように、本発明においては、危険率ｐ（ｘ，ｙ）をｘ，ｙの組み合わせそれぞれについて計算した後、重みｗ（ｐ（ｘ，ｙ））を求めるといった手順を踏む必要は必ずしもなく、危険率ｐ（ｘ，ｙ）を具体的に求めずに、合成関数としての（ｗ°ｐ）を、直接計算する構成としてもよい。
【００８６】
以上述べたように、分散σ^２の推定をし（例えば、（７）式のＥ[σ^２]）、かつ、パラメータＣを適当に決める（例えば、Ｃ＝３）ことにより、（４）式を用いて、ある画素ｘにつき定義される集合Ｎ（ｘ）（上述したように、例えば画素ｘを中心とした３×３画素分のエリア等）に含まれるすべての画素ｙについて、具体的な重みｗ１（ｐ（ｘ，ｙ））を求めることができる。後は、上記（２）式におけるｗ（ｐ（ｘ，ｙ））に代えて、このｗ１（ｐ（ｘ，ｙ））を用いることにより、コヒーレントフィルタの具体的な数値演算を実施することが可能となる。そしてその結果、時間分解能は勿論のこと、空間分解能をも損なわずに、ノイズを強く抑制した画素値ｖ´（ｘ）＝（ｖ´_１（ｘ），ｖ´_２（ｘ），…，ｖ´_Ｋ（ｘ））（＝（３）式）、すなわちそのようなＫ枚の静止画像ないし動画像を、得ることができる。
【００８７】
このような画像処理を、概念的に把握しやすいよう図示したものが、図３（ａ）乃至図３（ｃ）である。すなわちまず、図３（ａ）においては、１，２，…，Ｋ枚ある静止画像において、ある画素ｘにつき、該画素ｘを中心とした３×３画素分の矩形状エリアＮ^３×３（ｘ）が想定されている。この矩形状エリアＮ^３×３（ｘ）の左角隅における画素を、ｙ_１とすれば、この画素ｙ_１は、画素値ｖ（ｙ_１）を有している。
【００８８】
そして、この画素値ｖ（ｙ_１）を構成するスカラー値ｖ_１（ｙ_１），ｖ_２（ｙ_１），…，ｖ_Ｋ（ｙ_１）と画素値ｖ（ｘ）におけるスカラー値ｖ_１（ｘ），ｖ_２（ｘ），…，ｖ_Ｋ（ｘ）とのそれぞれにより、上記（４）式によって重みｗ１（ｐ（ｘ，ｙ_１））が計算される（図３（ｂ））。また、矩形状エリアＮ^３×３（ｘ）の残る画素ｙ_２，…，ｙ_８についても同様で、結局図３（ｂ）に示すように、ｗ１（ｐ（ｘ，ｙ_１）），…，ｗ１（ｐ（ｘ，ｙ_８））及び、ｗ１（ｐ（ｘ，ｘ））が得られる。この場合、（８）式より危険率ｐ（ｘ，ｘ）は、“１”であり、したがって重みｗ１（ｐ（ｘ，ｘ））も、（９）式より“１”である（＝最大の重み付けがされている）。
【００８９】
次に、このようにして得られた、重みｗ１（ｐ（ｘ，ｙ_１）），…，ｗ１（ｐ（ｘ，ｙ_８）），ｗ１（ｐ（ｘ，ｘ））を、対応する画素の、ｋ枚目の画像におけるスカラー値ｖ_ｋ（ｙ_１），ｖ_ｋ（ｙ_２），…，ｖ_ｋ（ｙ_８），ｖ_ｋ（ｘ）にそれぞれ乗算して総和を取り（上記（２）式における分子に該当する。）、これを矩形状エリアＮ^３×３（ｘ）に関する重みｗ１の総和（同じく（２）式の分母に該当する。）により除せば、当該ｋ枚目の画像における画素ｘについての、ノイズが抑制されたスカラー値ｖ´_ｋ（ｘ）を求めることができる（図３（ｃ））。また、ｋ＝１，２，…，Ｋのすべての画像につき、同じ重みｗ１（ｐ（ｘ，ｙ_１）），…，ｗ１（ｐ（ｘ，ｙ_８）），ｗ１（ｐ（ｘ，ｘ））を用いて、ノイズが抑制されたスカラー値ｖ´_ｋ（ｘ）を求めることによって、画素ｘにおけるノイズが抑制された画素値ｖ´_ｋ（ｘ）＝（ｖ´_１（ｘ），ｖ´_２（ｘ），…，ｖ´_Ｋ（ｘ））が得られる。すべての画素ｘにつき、上記演算を繰り返せば、ノイズを抑制したＫ枚の画像が得られる。
【００９０】
このようにしてコヒーレントフィルタで算出された画素値ｖ´（ｘ）で構成される画像では、オリジナル画像で見られたランダムなノイズが、十分に抑制される。
【００９１】
なお、以上までに述べた各処理は、例えば図４（ａ），図４（ｂ）に示すようなフローチャートに則ってこれを行えばよく、また、当該各処理に係る演算・画像表示等を実際のＸ線ＣＴ装置１００上で実現するためには、例えば、図２に示すように、分散値推定部１１１、重み演算部１１２及び画素値演算部１１３により構成される画像処理部１１０を設けて、これを実施すればよい。
【００９２】
このうち重み演算部１１２は、上述した手順通り、画素値ｖ（ｘ）及びｖ（ｙ）から直接重みｗ１（ｐ（ｘ，ｙ））を求める構成となっている。したがって当該演算部１１２は、危険率ｐ１（ｘ，ｙ）の値を具体的に求めることなく、重みを直接に求める装置である。なお、上記したような構成ではなく、具体的に危険率ｐ１（ｘ，ｙ）の値を求める危険率演算部（適合度定量化部）と、その出力に基づいて重みｗ１（ｐ（ｘ，ｙ））を求める重み演算部という、二段の手順を踏む構成としてもよい。いずれにせよ、重み演算部１１２は、分散値推定部１１１により推定された分散σ^２と、ｖ（ｘ）及びｖ（ｙ）を用いて重みｗ１（ｐ（ｘ，ｙ））を算出する。
【００９３】
また、画素値演算部１１３は、画素値ｖ（ｘ）及びｖ（ｙ）、並びに重み演算部１１２により数値演算された重みｗ１（ｐ（ｘ，ｙ））を使って、画素値ｖ´（ｘ）を演算する。すなわち当該演算部１１３は、元となる画像のノイズを抑制する処理、すなわちコヒーレントフィルタの適用を実際に行う（以下、これを「コヒーレントフィルタをかける」と表現する）。
【００９４】
上記のようなダイナミック・コヒーレントフィルタ処理においてＫ枚の静止画像から構成される動画像に、コヒーレントフィルタをかける場合には、上記画像処理部１１０における処理は、一旦すべての静止画像を再構成した後、これらを上記記憶装置１０Ｍに蓄え、後処理として後にこれらに対してコヒーレントフィルタをかけるようにしてもよいが、本実施形態はこのような形態に限定されるものではなく、上述した連続スキャン、連続投影データ収集、連続再構成及び連続表示という流れの中で、コヒーレントフィルタをかける処理をリアルタイムに実施する（以下、これを「リアルタイム・コヒーレントフィルタ処理」と呼ぶ。）のでもよい。
【００９５】
リアルタイム・コヒーレントフィルタ処理の好ましい実施形態においては、新しい画像が撮影され再構成されるたびに、以下のような処理を行う。最初に得られた画像（画像番号１）から最新の画像（画像番号Ｍ）までのうち、画像番号Ｍ，Ｍ−１，…，Ｍ−Ｋ＋１を持つＫ枚の静止画像上、共通の同一点（同一座標）の画素ｘの持つ画素値（スカラー値）を並べてＫ次元ベクトル値ｖ（ｘ）＝（ｖ_Ｍ（ｘ），ｖ_Ｍ−１（ｘ），…，ｖ_{Ｍ−Ｋ＋１}（ｘ））を構成する。こうして、上記の「ダイナミック・コヒーレントフィルタ処理」と全く同様にコヒーレントフィルタをかけることができる。ただし、画素値演算部１１３は実際には画素値ｖ´（ｘ）の全ての要素を計算するのではなく、最新の画像（画像番号Ｍ）に対応するスカラー値ｖ_Ｍ´（ｘ）だけを計算する。この結果、計算速度が向上するので、リアルタイムでノイズが抑制された最新の画像を表示できる。
【００９６】
この「リアルタイム・コヒーレントフィルタ処理」の別の好ましい実施形態として、最初のＫ枚の画像が得られた時点で、上記と全く同様にコヒーレントフィルタをかけてｖ_１´（ｘ），…，ｖ_Ｋ´（ｘ）を求めておき、以後は、Ｋ次元ベクトル値を画像番号Ｍ，Ｍ−１，…，Ｍ−Ｋ＋１を持つＫ枚の静止画像を用いてｖ（ｘ）＝（ｖ_Ｍ（ｘ），ｖ_Ｍ−１´（ｘ），…，ｖ_{Ｍ−Ｋ＋１}´（ｘ））によって構成し、これに対して上記のリアルタイム・コヒーレントフィルタ処理を適用するように構成してもよい。なお、これらのリアルタイム・コヒーレントフィルタ処理の際に画素値ベクトルｖ（ｘ）の次元Ｋを、マニュアル設定、あるいは自動設定によって、随時変更できるように構成しておくと便利である。
【００９７】
このようにコヒーレントフィルタにより、空間及び時間分解能を低下させることなく、ノイズだけを効果的に抑制したＣＴ画像を使ってＣＢＰスタディを実行し、組織内の局所的な血流動態、つまり局所組織内の毛細血管を通過する血流の動態を定量的に解析し、その局所血流動態を表すインデックス（ＣＢＰ、ＣＢＶ、ＭＴＴ、Ｅｒｒ）を求めることにより、その精度及び信頼性の向上が期待できる。
【００９８】
以上のように分解能の低下を抑え、ノイズを除去した画像に対してＣＢＰスタディ処理が実施される。
（ＣＢＰスタディ）
上述したように、ＣＢＰスタディでは、脳組織内の”毛細血管を通過する血流”の動態を定量的に表すＣＢＰ、ＣＢＶ、ＭＴＴ、Ｅｒｒのインデックスを求め、またこれらインデックスの空間的分布を表すマップを出力する。
【００９９】
ＣＢＰ：脳組織の毛細血管内の単位体積及び単位時間あたりの血流量［ml/１００ml/miｎ］
ＣＢＶ：脳組織内の単位体積あたりの血液量［ml/１００ml］
ＭＴＴ：毛細血管の血液平均通過時間［秒］
Ｅｒｒ：解析モデルからの実測値のずれ残差の指標
なお、残差Ｅｒｒが低いことは、参照脳動脈から支配を受けている可能性が高いことを意味し、逆に、残差Ｅｒｒが高いことは、参照脳動脈から支配を受けている可能性が低いことを意味している。
【０１００】
ＣＢＰスタディでは、トレーサーとして脳血管透過性を持たない造影剤、たとえばヨード造影剤が使用される。インジェクターにより肘静脈から急速に注入されたヨード造影剤は、心臓、肺を経由して、脳動脈から流れ込む。そして、脳動脈から、脳組織内の毛細血管を経て、脳静脈へと流れ出ていく。このとき、脳血管透過性を持たない造影剤、たとえばヨード造影剤は正常な脳組織内の毛細血管では造影剤は血管外へ漏れ出ることなく通過する。
【０１０１】
造影剤の通過の様子をダイナミックＣＴで連続的に撮影して、その連続画像から、脳動脈上の画素の時間濃度曲線Ｃａ（ｔ）、毛細血管を含む脳組織上の画素の時間濃度曲線Ｃｉ（ｔ）、脳静脈上の画素の時間濃度曲線Ｃｓｓｓ（ｔ）をそれぞれ測定する。
【０１０２】
ＣＢＰスタディでは、造影剤の血中濃度について脳組織に近い脳血管の血中濃度の時間曲線Ｃａ（ｔ）と、毛細血管の血中濃度の時間曲線Ｃｉ（ｔ）との間で成り立つ理想的な関係を解析モデルとして採用する。つまり脳組織に入る直前の血管から造影剤を注入した場合、毛細血管を含む脳組織単位体積（１画素）内の時間濃度曲線は垂直に立ち上がり、一定値を維持し、そして若干の勾配を持って立ち下がる。これは、矩形関数で近似することができる（ｂｏｘ−ＭＴＦ法：ｂｏｘ−ＭｏｄｕｌａｔｉｏｎＴｒａｎｓｆｅｒＦｕｎｃｔｉｏｎｍｅｔｈｏｄ）。
【０１０３】
脳動脈血中時間濃度曲線Ｃａ（ｔ）を入力関数、脳組織の時間濃度曲線Ｃｉ（ｔ）を出力関数として、毛細血管を通過する過程の特徴を、矩形関数で表される伝達関数として求めることができる。
【０１０４】
（具体的な手順）
図５、図６には、本実施形態によるＣＢＰスタディの典型的手順を示している。まず、肘静脈等の血管にボーラスインジェクション（造影剤を一気に投与する）を行い、その直後あるいは直前からダイナミックＣＴ（同じ箇所を反復して撮影する）を行う。最も典型的な手技として、肘静脈へボーラスインジェクションを行った場合、概ね２０〜４０秒の間、例えば０．５〜２秒間隔で撮影を繰り返す。ダイナミックＣＴで得たＮ枚のＣＴ画像のうちのｊ枚目の各ピクセル（ｘ、ｙ）のＣＴ値をｖ（ｘ、ｙ、ｊ）とする。これはこの画素（ｘ、ｙ）における時間濃度曲線（滑らかな曲線である）ｆ（ｔ、ｘ、ｙ）をサンプリングしたものに他ならない。
【０１０５】
まず、前処理として、ステップＳ１で、ＣＴ画像各々から、明らかに脳組織以外の組織であることが判別される画素を、解析対象から除外する。すなわち、脳組織のＣＴ値として考えられる範囲（例えばＣＴ値１０〜６０ＨＵ）に入らない値を示す画素は、空気や骨、脂肪などに対応する画素であり、脳血流の定量とは関係ないのでこれらは無視して良い。この解析範囲は、デフォルトとして、１０〜６０ＨＵに設定されるが、入力部１０９を介して任意に設定可能である。
【０１０６】
また、前処理として、ステップＳ２で、造影効果の初期化が行われる。各画素に於ける造影効果（ＣＴ値の上昇）を得るためには、各画素（ｘ，ｙ）について、その画素に対応する組織に造影剤が到達する以前の画像（一般に複数枚得られる）を、通し番号１，２，…Ｋで表すと、その時間的平均値は、
【数８】

を求め、この値をｂ（ｘ，ｙ）とする。そして、ｊ＝Ｋ＋１、Ｋ＋２、…，Ｎの各画像の画素値ｖ（ｘ，ｙ，ｊ）について、
ｑ（ｘ，ｙ，ｊ）＝ｖ（ｘ，ｙ，ｊ）−ｂ（ｘ，ｙ）
ｊ＜Ｋについて
ｑ（ｘ，ｙ，ｊ）＝０
とすればよい。処理を簡単にするためには、どの画素に関しても同じＫを採用しても良い。こうして得られたｑ（ｘ，ｙ，ｊ）は、滑らかな連続曲線の時間濃度曲線ｑ（ｔ，ｘ，ｙ）をｔ＝ｔ１，ｔ２，…ｔＮにおいてサンプリングしたものに他ならないと考えることができる。このｑ（ｔ，ｘ，ｙ）を用いて脳血流の定量解析を行う。
【０１０７】
定量解析にあたってはまず、右脳エリアと左脳エリアをＣＴ画像上で分離することができる。上述したようにＣＢＰスタディでは、毛細血管の血流動態の様子を、脳動脈の時間濃度曲線Ｃａ（ｔ）に対する脳組織の時間濃度曲線Ｃｉ（ｔ）の伝達関数ＭＴＦとして求めるものであり、従って、解析対象の脳組織が、参照曲線Ｃａ（ｔ）の脳動脈の支配下にないならば、計算しても無駄である。少なくとも左脳と右脳とでそれぞれ別々の脳動脈の時間濃度曲線Ｃａ（ｔ）を使って個別に解析する、つまり左脳の脳動脈の時間濃度曲線Ｃａ（ｔ）は同じ左脳の脳組織の解析にだけ使用し、同様に、右脳の脳動脈の時間濃度曲線Ｃａ（ｔ）は同じ右脳の脳組織の解析にだけ使用することは無駄な計算を減らす効果がある。
【０１０８】
脳を左脳エリアと右脳エリアとを分割するために、図７に示すように、ＣＴ画像上に分割線が画面上に図形として重ねて表示される（Ｓ３）。分割線が最初は画像中央に表示されるように構成しても良い。操作者は、画像を参照して、分割線を移動し、また分割線を構成する複数の構成点を移動して任意に屈曲させることにより、左右エリアを分割する。
【０１０９】
このように脳を左脳エリアと右脳エリアとに分けて、それぞれのエリアに解析範囲を限局することで、解析処理工数を減らすことができる。つまり左脳の脳動脈（左ＡＣＡ、左ＭＣＡ、左ＰＣＡ）の時間濃度曲線Ｃａ（ｔ）は左脳エリアの解析（伝達関数最適化処理）にだけ使用し、同様に、右脳の脳動脈（右ＡＣＡ、右ＭＣＡ、右ＰＣＡ）の時間濃度曲線Ｃａ（ｔ）は同じ右脳の脳組織の解析にだけ使用する。解析処理工数を減らすために、左脳エリアと右脳エリアをそれぞれさらに領域分割し、さらに狭いエリアに解析処理を限局しても良い。
【０１１０】
この領域分割には、幾何学図法、例えばボロノイ図法が採用される。ボロノイ図法は、周知の通り、病院、店舗、消防署等の施設の最適な配置等の分野によりよく用いられる手法であり、平面上に配置された多数点（店舗等に相当、母点）からの距離に応じて平面を複数の勢力域に分割することを特徴とする。
【０１１１】
図８に示すように、本実施形態では、ボロノイ図法は、左脳エリアと右脳エリアとで個別に適用される。左脳エリアは、左ＡＣＡ、左ＭＣＡ、左ＰＣＡを３つの母点として、左ＡＣＡの勢力域と、左ＭＣＡの勢力域と、左ＰＣＡの勢力域とに分割される。左ＡＣＡ、左ＭＣＡ、左ＰＣＡに対応する３つの母点を通過する円の中心にボロノイ点を設定する。ボロノイ点を中心として、左ＡＣＡと左ＭＣＡの２つの母点の垂直二等分線と、左ＭＣＡと左ＰＣＡの２つの母点の垂直二等分線と、左ＡＣＡと左ＰＣＡの２つの母点の垂直二等分線とが連結される。これら垂直二等分線により左脳エリアが３つの勢力域に分割される。同様に、右脳エリアは、右ＡＣＡ、右ＭＣＡ、右ＰＣＡを３つの母点として、右ＡＣＡの勢力域と、右ＭＣＡの勢力域と、右ＰＣＡの勢力域とに分割される。
【０１１２】
左ＡＣＡの時間濃度曲線Ｃａ（ｔ）に対する脳組織の時間濃度曲線Ｃｉ（ｔ）の伝達関数ＭＴＦを、左ＡＣＡの勢力域に限局して画素ごとに求める。同様に、左ＭＣＡ、左ＰＣＡ、右ＡＣＡ、右ＭＣＡ、右ＰＣＡに対してそれぞれの勢力域に限局して画素ごとに伝達関数ＭＴＦを求める。
【０１１３】
このように左脳エリアと右脳エリアとをそれぞれ複数の勢力域に分割して、それぞれの勢力域に解析範囲を限局することで、解析処理工数をさらに減らすことができる。
【０１１４】
次に、ＣＴ画像上で脳動脈上に脳動脈ＲＯＩを設定するのであるが、この設定の精度を向上し且つ容易にするためにＲＯＩ設定支援部１２１により支援マップが作成され、ＣＴ画像とは別に、またはＣＴ画像に重ねて表示される（Ｓ４）。支援マップとしては、図９（ａ）乃至図９（ｃ）に示すように、例えば、ＡＴ（アピアランスタイム）マップ、ＰＴ（ピークタイム）マップ、ＴＴ（トランジットタイム）マップがあげられる。各画素について、図１０に示すように、造影前の任意の時刻（例えばデータ収集開始時刻）Ｔ０から造影剤濃度がピークｐｅａｋの数パーセント（例えば１パーセント）に達する時刻までの時間ＡＴ、時刻Ｔ０から造影剤濃度がピークに達した時刻までの時間（ピークタイム）ＰＴ、又は造影剤の移動時間を例えば半値幅で表すＴＴが計算され、マップとして生成され表示される。デフォルトでは、これらＡＴマップ、ＰＴマップ、ＴＴマップの全種類が生成され表示されるようになっているが、任意の１種類、又は任意の２種類を操作者が選択することが可能である。
【０１１５】
これら数値は脳動脈では他の組織と比較して高い値で現れる傾向にあるので、その値を中心とした値を持つ画素だけを表示するよう設定されたカラールックアップテーブルを通してカラー表示させることで、脳動脈の場所を容易に識別し、脳動脈ＲＯＩを正確に設定することが可能となる（Ｓ５）。典型的には、左右脳エリアそれぞれに、脳動脈ＲＯＩは、前大脳動脈（ＡＣＡ）、中大脳動脈（ＭＣＡ）、後大脳動脈（ＰＣＡ）の３箇所ずつ設定される。
【０１１６】
なお、マルチスライス、撮影の場合、例えば隣接する４枚のスライスを解析対象とする場合、図１１に示すように、その各スライスで個々に脳動脈ＲＯＩを設定することは作業負担が大きいばかりで、解析を行う上では必要のない作業である。従って、ある任意の１スライスで設定した脳動脈ＲＯＩを他のスライスにも共用する。または、後述するコヒーレントレグレッション法を使って全スライスで共通に用いることのできる脳動脈の時間濃度曲線Ｃａ（ｔ）を作成するようにしてもよい。
【０１１７】
次に、設定された脳動脈ＲＯＩ各々に関して時間濃度曲線Ｃａ（ｔ）が、ダイナミックＣＴによる連続画像データから時間濃度曲線作成部１２２により作成される（Ｓ６）。
【０１１８】
ここで、画素サイズに比べて、脳動脈は非常に細いものが多く、しかも一般にＣＴの撮影スライスに対して直交していないため、画像上のどの画素も正確に動脈血のＣＴ値を表しておらず、１画素が脳動脈と他の組織との混在により構成され、そのパーシャルボリューム効果のためにそれよりも低い造影効果しか示さないことがほとんどである。また、動脈をパーシャルボリュームとして含むこれらの画素において、画像ノイズが大きい。特に脳梗塞を生じている部位の動脈等においては、造影効果が比較的小さいために、ノイズの影響は甚大である。画像ノイズに関しては上述したコヒーレントフィルタにより抑制されているものの、パーシャルボリューム効果の影響は依然として残存している。
【０１１９】
この問題は、脳動脈の時間濃度曲線は、単一のスライス画像において計測するのではなく、その動脈を含む立体内の画素を用いて後述するコヒーレントレグレッション法を適用することで抑制することが可能である。従って、上述したコヒーレントフィルタ法に代えて、この段階でコヒーレントレグレッション法を適用するようにしてもよい。
【０１２０】
また、この方式によれば、動脈ごとにそれに対応する唯一のスライス像の時間濃度曲線が得られ、従って、撮影範囲内にある全スライス中の任意の部位の解析に利用することができ、これによって、特定の動脈について、その脳動脈の時間濃度曲線が最も明瞭に得られるスライスを選んで、その脳動脈の時間濃度曲線を全スライスに適用することができ、脳動脈の時間濃度曲線の数を減らすことができる。
【０１２１】
（コヒーレントレグレッション法）
上記時間濃度曲線の作成では、パーシャルボリューム効果、ランダムノイズの影響を除去することが重要である。まず、「時間濃度曲線」とは、上記ダイナミックＣＴ画像中の特定の部位におけるＣＴ値（濃度値）の経時的変化を表す曲線である。ことに、上記医用画像診断装置においては、人体組織等における血流動態や代謝機能等の詳細を調べる事を目的として、人体の特定組織内の造影剤濃度等の経時的変化を時間濃度曲線として計測することが行われている。また、天体観測等においては、特定の天体の光度変化等を解析する目的で、時間濃度曲線が用いられる。より形式的に明示すると、すなわち、時間濃度曲線とは、時刻ｔ_ｋにおけるある部位の濃度値をｄ_ｋとするとき、対の列｛＜ｔ_ｋ，ｄ_ｋ＞（ｋ＝１，２，・・・，Ｋ）｝として表現される。また、時間濃度曲線の多くの用途においては、必ずしもｄ_ｋの絶対的な値が必要なのではなく、むしろ最初の画像１を基準とする増分（ｄ_ｋ−ｄ_１）だけが得られれば十分である。さらにそのような用途のうちの多くでは、単に（ｄ_ｋ−ｄ_１）に比例するデータＡ（ｄ_ｋ−ｄ_１）（ここにＡは未知の定数）だけが得られれば十分である。この場合には、従って、対の列｛＜ｔ_ｋ，Ａ（ｄ_ｋ−ｄ_１）＞（ｋ＝１，２，・・・，Ｋ）｝が、求める時間濃度曲線である。
【０１２２】
このような時間濃度曲線を求めるためには、原理的には、上記ダイナミックＣＴ画像を構成する各画像ｋ（ｋ＝１，２，・・・，Ｋ）における、該時間濃度曲線を測定しようとする部位に含まれる画素ｘのスカラー値ｖ_ｋ（ｘ）を用いて、対の列｛＜ｔ_ｋ，ｖ_ｋ（ｘ）＞｝あるいは、｛＜ｔ_ｋ，Ａ（ｖ_ｋ（ｘ）−ｖ_１（ｘ））＞｝を構成すればよい。
【０１２３】
しかし、実用においては、上記医用画像診断装置等によって撮影されたダイナミックＣＴ画像にランダムなノイズが含まれているために、本来測定しようとする時間濃度曲線を正確に求められないという問題がある。
【０１２４】
さらに、実用においては、これらのダイナミックＣＴ画像においては、いわゆる「パーシャルボリューム効果」が生じる。パーシャルボリューム効果とは、すなわち、被検体内の微小な物体の像は、画像上では少数個の画素によって表現されるが、これら少数個の画素には、被検体内の隣接する物体の像も影響を与えるため、これら少数個の画素の画素値は（本来計測しようとする濃度値の変動に比例するものの）比較的小さな変動しか示さない、という現象である。言い換えれば、これら少数個の画素の画素値は僅かな信号しか含まない。従って、パーシャルボリューム効果が生じている場合には、どの画素ｘを取っても対の列｛＜ｔ_ｋ，ｖ_ｋ（ｘ）＞（ｋ＝１，２，・・・，Ｋ）｝は非常に信号レベルが低く、本来計測しようとしているのではない組織における濃度値の変化の影響を受け、さらにランダムなノイズが存在するために、本来測定しようとする時間濃度曲線｛＜ｔ_ｋ，ｄ_ｋ＞｝を正確に求められないという問題がある。
【０１２５】
そこで、従来は、ランダムなノイズを抑制するために、時間的又は空間的平滑化を用いていたが、時間平均を行うと時間分解能が損なわれてしまい、また、空間平均を行うと、本来の該測定しようとする部位以外の部位の濃度の経時変化が計測値に混入するという問題点があった。このような問題点を解決し、より正確な時間濃度曲線を得るために、コヒーレントフィルタを採用する。
【０１２６】
まず、本実施形態のコヒーレントフィルタにおいて用いるべき、帰無仮説について説明する。計測しようとする部位における真の時間濃度曲線を｛＜ｔ_ｋ，ｄ_ｋ＞（ｋ＝１，２，・・・，Ｋ）｝であると仮定するとき、その一次変換である｛＜ｔ_ｋ，Ａ（ｄ_ｋ−ｄ_１）＞（ｋ＝１，２，・・・，Ｋ）｝（だたしＡは未知の係数）を計測することを目的とする場合において、計測しようとする部位に概ね相当する画素の集合Ｒを設定する。この集合Ｒの要素である任意の画素ｘ∈Ｒについて、条件Ｑ：「もし、この画素ｘが上記真の時間濃度曲線を良く反映し、しかも他の部位の経時的濃度変化の影響をほとんど受けていない」のであれば、ベクトル値としての画素値ｖ（ｘ）＝（ｖ_１（ｘ），ｖ_２（ｘ），．．．，ｖ_Ｋ（ｘ））について、パーシャルボリューム効果およびランダムノイズの影響を考慮することによって、
ｖ_ｋ（ｘ）＝ｐ（ｘ）ｄ_ｋ＋ｑ（ｘ）＋γ_ｋ（ｘ） …（１１）
（ｋ＝１，２，・・・，Ｋ）
が成り立つと仮定することができる。ここに、ｐ（ｘ）およびｑ（ｘ）は、画素ｘごとに異なるが画像番号ｋ（すなわち撮影時刻ｔ_ｋに相当）によっては変化しない未知の係数であり、パーシャルボリューム効果をモデル化したものである。またγ_ｋ（ｘ）はランダムなノイズをモデル化したものであって、画素ｘごとに、しかも画像番号ｋごとに値が異なるが、その期待値は０であり、またその統計分布は画素ｘにも画像番号ｋにも依存しない。
【０１２７】
以上の仮定によれば、該集合Ｒの要素である任意の２個の画素ｘ，ｙに関して、もし「画素ｘ，ｙが共に上記の条件Ｑを満たす。」という命題が成り立つのであれば、次式の関係が成り立つことが証明できる。
【０１２８】
ｖ_ｋ（ｘ）＝ａ_１ｖ_ｋ（ｙ）＋ａ_２＋ξ_ｋ（ｋ＝１，２，・・・，Ｋ）
… （１２）
ここに、ａ_１およびａ_２は、画素の組ｘ，ｙごとに異なるが画像番号ｋ（すなわち撮影時刻ｔ_ｋ）によっては変化しない未知の係数である。またξ_ｋはランダムなノイズであって、画素の組ｘ，ｙごとに、しかも画像番号ｋごとに値が異なるが、その期待値は０である。
【０１２９】
（１２）式は以下のようにして導かれる。すなわち、ｘにｙを代入して得られる式
ｖ_ｋ（ｙ）＝ｐ（ｙ）ｄ_ｋ＋ｑ（ｙ）＋γ_ｋ（ｙ）… （１３）
を変形すると、
【数９】

とおくことによって、（１２）式が導かれる。ここで、（１６）式のａ_１とａ_２は、パーシャルボリューム効果を現すパラメータであり、また（１６）式のξ_ｋは、ランダムなノイズを表す。
【０１３０】
以上から、「画素ｘ，ｙが共に条件Ｑを満たす。」という命題は、帰無仮説Ｈ_０´「ｖ_ｋ（ｘ）＝ａ_１ｖ_ｋ（ｙ）＋ａ_２＋ξ_ｋ（ｋ＝１，…，Ｋ）である。」と等価であることが示された。
【０１３１】
次に帰無仮説Ｈ_０´「ｖ_ｋ（ｘ）＝ａ_１ｖ_ｋ（ｙ）＋ａ_２＋ξ_ｋ（ｋ＝１，…，Ｋ）である。」を、実質的に等価であり、かつ実際に検定できる形式の命題に変換する方法について述べる。この帰無仮説を改めて数学的に厳密な表現で述べると、帰無仮説Ｈ_０´「ある定数ａ_１およびａ_２が存在して、ξ_ｋ＝ｖ_ｋ（ｘ）−ａ_１ｖ_ｋ（ｙ）−ａ_２（ｋ＝１，…，Ｋ）は平均０、分散（σ^２ｈ（ａ_１））の正規分布に従う。」となる。ここに係数ｈ（ａ_１）は、
ｈ（ａ_１）＝１＋ａ_１ ^２ … （１７）
である。（１７）式はａ_１とξ_ｋの定義である（１６）式、および、ランダム変数に関する分散の持つ一般的な性質から直ちに導かれる。また、上記の分散σ^２の値は、簡便かつ実用上十分正確に推定できる。
【０１３２】
以上から、もし、上記の定数ａ_１およびａ_２を決定することができれば、上記の帰無仮説Ｈ_０´を検定することが可能である。そして実際上は、これらの定数の最適な推定値ａ_１ ^〜およびａ_２ ^〜が得られれば十分である。
【０１３３】
このような、定数ａ_１およびａ_２の最適な推定値の算出には、公知の当てはめ法（fitting）がそのまま利用できる。そこで、以下では、そのような当てはめ法の典型的な具体例として、線形最小二乗法を用いる場合における概要を説明する。線形最小二乗法を本実施形態に適用するには、単に、上記の帰無仮説のξ_ｋの二乗和をＳ（ａ）として、すなわち
【数１０】

を定義する。Ｓ（ａ）の値は定数ベクトルａ＝（ａ_１，ａ_２）、すなわち上記の定数ａ_１およびａ_２、の値に依存する。このＳ（ａ）が最小の値を取るような定数ベクトルａを算出すれば、定数ａ_１およびａ_２に関する、不偏推定の意味での最適な推定値ａ_１ ^〜およびａ_２ ^〜が得られる。なお、線形最小二乗法の具体的な計算方法としては、様々な公知の方法を利用することができ、しかも、これら公知の計算方法はいずれも非常に簡単であり、必要な計算時間はごく僅かである。
【０１３４】
このようにして、上記の定数ａ_１，ａ_２の最適な推定値ａ_１ ^〜，ａ_２ ^〜を算出した結果、次式で定義される残差
ｒ_ｋ ^〜（ｘ，ｙ）＝ｖ_ｋ（ｘ）−ａ_１ ^〜ｖ_ｋ（ｙ）−ａ_２ ^〜 …（１９）を具体的に計算することができる。従って、この残差ｒ_ｋ ^〜を用いて、上記の帰無仮説Ｈ_０´を、実質的に等価な帰無仮説Ｈ_０”「ｒ_ｋ ^〜（ｘ，ｙ）（ｋ＝１，…，Ｋ）は平均０、分散（１＋（ａ_１ ^〜）^２）σ^２の正規分布に従う。」と言い換えることができる。これは、実際に検定の計算を実行可能な具体的命題である。
【０１３５】
なお、さらに、ベクトルによる表現
【数１１】

（ただし、ベクトルａ及びξは画素の組ｘ，ｙに依存する。）を導入し、また、次式
ｆ（ａ^~，ｖ（ｙ））＝ａ_１ ^〜ｖ（ｙ）＋ａ_２ ^〜 … （２１）
で定義されるベクトル値関数ｆを用いて帰無仮説Ｈ_０´を言い換えると、帰無仮説Ｈ_０”は「ｖ（ｘ）＝ｆ（ａ^~，ｖ（ｙ））＋ξ、（ただし、ξは平均０、分散（１＋（ａ_１ ^〜）^２）σ^２の正規分布に従う。）」となり、これは上述した帰無仮説Ｈ_０と全く同じ形式である。すなわち、本実施形態は上述したコヒーレントフィルタの一変形例であることは明らかである。なお、ここで、上記ｆ（ａ^~，ｖ（ｙ））とは、すなわち、画素ｙの画素値ｖ（ｙ）に対して、パーシャルボリューム効果を現すパラメータａを最適に調節して、画素ｘの画素値ｖ（ｘ）と最も高い適合度を持つように変換したものを意味する。
【０１３６】
次に、本実施形態において、上記の帰無仮説Ｈ_０”を用いて、コヒーレントフィルタによって時間濃度曲線を求める方法について説明する。計測しようとする部位に概ね相当する画素の集合Ｒについて、この集合Ｒに含まれるあるひとつの画素ｘ∈Ｒについて、集合Ｒの要素である全ての画素y∈Ｒに対して、以下の計算を行う。すなわち、上記の方法を用いて実際に残差ｒ_ｋ ^〜（ｘ，ｙ）（ｋ＝１，…，Ｋ）を算出し、次に、上記の帰無仮説Ｈ_０”「ｒ_ｋ ^〜（ｘ，ｙ）（ｋ＝１，…，Ｋ）は平均０、分散（１＋（ａ_１ ^〜）^２）σ^２の正規分布に従う。」を棄却する場合の危険率ｐ（ｘ，ｙ）ないし重みｗ（ｐ（ｘ，ｙ））を具体的に計算する。そして、重み付き平均ｖ´_ｋ（ｘ）を下式（２２）によって計算し、画素ｘにおける時間濃度曲線｛＜ｔ_ｋ，ｖ´_ｋ（ｘ）−ｖ´_１（ｘ）＞（ｋ＝１，２，・・・，Ｋ）｝を構成する。
【０１３７】
【数１２】

【０１３８】
こうして得られた時間濃度曲線は、画素ｘにおける真の時間濃度曲線｛＜ｔ_ｋ，ｄ_ｋ＞｝の一次変換である｛＜ｔ_ｋ，Ａ（ｄ_ｋ−ｄ_１）＞｝（だたしＡは未知の係数）を近似している計測値であり、しかも、重み付き平均の効果によって、ランダムなノイズが抑制されている。また、他の画素ｙの画素値ベクトルに対しては、式から明らかなように、パーシャルボリューム効果の影響を補正したものが用いられている。さらに、本実施形態はコヒーレントフィルタの共通の特徴である「時間平均を全く使用せず、また空間平均を画素ｘとの適合度に基づく重みを使って計算する」という性質を有する。従って、本実施形態によって、時間分解能を損なわず、パーシャルボリューム効果の影響を抑制し、しかもランダムなノイズが抑制された時間濃度曲線を得ることができる。なお、このようにして時間濃度曲線を求める方式を、特に「コヒーレントレグレッション法」と称す。
【０１３９】
次に、具体的に、医療用のＸ線ＣＴにおけるダイナミックＣＴ撮影等で得られたダイナミックＣＴ画像における、時間濃度曲線の臨床的利用の一例を説明する。この応用例では、造影剤を血管に急速に注入しながら、ダイナミックＣＴ等の撮影を行い、人体組織中に存在する動脈の像の濃度変化を時間濃度曲線として計測することによって、当該組織における血流動態を診断しようとするものである。
【０１４０】
この応用例において、多くの場合、人体組織中の動脈は一般に非常に細いために、ＣＴによる断層画像上に現れる動脈の像は、パーシャルボリューム効果を生じる。さらに、像にはランダムなノイズが含まれていることは言うまでもない。このため、従来の方法では、動脈に関する十分に正確な時間濃度曲線を得ることは困難であり、強いて計測を行えば、動脈に関する真の時間濃度曲線＜ｔ_ｋ，Ｄ_ｋ＞の一次変換である＜ｔ_ｋ，Ａ（Ｄ_ｋ−Ｄ_１）＞（ここにＤ_ｋは動脈の像に相当する一群の画素の、時刻ｔ_ｋにおける（スカラー値である）画素値を表す。また、ｋ＝１，２，・・・，Ｋ）をある程度近似する測定値＜ｔ_ｋ，（ｖ_ｋ（ｘ）−ｖ_１（ｘ））＞しか得られなかった。この測定値はランダムなノイズを含む。また、パーシャルボリューム効果の影響のために、係数Ａは未知のままである。
【０１４１】
そこで、＜ｔ_ｋ，Ａ（D_ｋ−D_１）＞を十分に近似する測定値＜ｔ_ｋ，（ｖ´_ｋ（ｘ）−ｖ´_１（ｘ））＞（ｋ＝１，２，・・・，Ｋ）を得ることができる。一方、同じ断層画像上で観察できる静脈の中には、相当に太いものが存在し、従ってそれらの静脈に関しては、従来の方法で、時間濃度曲線の十分に良い近似値＜ｔ_ｋ，（Ｊ_ｋ−Ｊ_１）＞（ｋ＝１，２，・・・，Ｋ）を得ることができる。ここにＪ_ｋは静脈の像に相当する一群の画素の、時刻ｔ_ｋにおける画素値を表す。
【０１４２】
ところで、血液循環に関する時間濃度曲線においては、命題Ｓ:「もし、時刻ｔ_１における血中の造影剤濃度が０であるならば、どの血管ｄに関する時間濃度曲線＜ｔ_ｋ，（ｄ_ｋ-ｄ_１）＞も、その曲線下面積（ＡＵＣ：Area Under Curve）が一致する」という性質が成り立つことが知られている。ここで言う曲線下面積とは、時間濃度曲線＜ｔ_ｋ，（ｄ_ｋ-ｄ_１）＞の時間tに関する積分を意味する。
【０１４３】
従って、ある血管ｄに関する時間濃度曲線＜ｔ_ｋ，（ｄ_ｋ-ｄ_１）＞の曲線下面積ＡＵＣ（ｄ）は、例えば次式によって近似的に計算することができる。
【０１４４】
【数１３】

【０１４５】
従って、静脈に関して従来の方法で得られた時間濃度曲線｛＜ｔ_ｋ，（Ｊ_ｋ-Ｊ_１）＞｝に関する曲線下面積ＡＵＣ（Ｊ）を（２２）式を用いて計算することができる。（ｄにＪを代入すればよい。）また、動脈に関して、仮に、時間濃度曲線｛＜ｔ_ｋ，（D_ｋ−D_１）＞｝が知られていれば、曲線下面積ＡＵＣ（D）を（１８）式を用いて同様に計算することができ、しかも上記命題Ｓに従って、
ＡＵＣ（D）≒ＡＵＣ（Ｊ） … （２４）
が成り立つ筈である。しかし、実際には、時間濃度曲線＜ｔ_ｋ，（D_ｋ−D_１）＞は未知であるため、ＡＵＣ（D）は計算できない。
【０１４６】
一方、本実施形態に係る方式で得られた時間濃度曲線＜ｔ_ｋ，（ｖ´_ｋ（ｘ）−ｖ´_１（ｘ））＞は、＜ｔ_ｋ，Ａ（D_ｋ−D_１）＞を近似するものであり、後者は未知の係数Ａを含んでいる。このため、｛＜ｔ_ｋ，（ｖ´_ｋ（ｘ）−ｖ´_１（ｘ））＞｝から（２３）式を用いて具体的に計算できる曲線下面積ＡＵＣ（ｖ´）は、ＡＵＣ（D）のちょうどＡ倍でなくてはならない。すなわち、
ＡＵＣ（ｖ´）≒Ａ・ＡＵＣ（D） … （２５）
である。すなわち、（２４）式と（２５）式から、
Ａ≒ＡＵＣ（ｖ´）／ＡＵＣ（Ｊ） … （２６）
という関係が成り立つ。（２６）式の右辺は（２３）式を用いて具体的に計算できるため、未知であった係数Ａの値が具体的に決定できる。そこで、この係数Ａの値を用いて時間濃度曲線＜ｔ_ｋ，（ｖ´_ｋ（ｘ）−ｖ´_１（ｘ））／Ａ＞を構成すれば、これは、動脈の時間濃度曲線＜ｔ_ｋ，（D_ｋ−D_１）＞を近似するものに他ならない。このように、曲線下面積を用いて、未知であった定数Ａの値を決定した時間濃度曲線を構成する方法を「ＡＵＣ法」と呼ぶ。
【０１４７】
以上から、ダイナミックＣＴ撮影等で得られたダイナミックＣＴ画像における、時間濃度曲線の臨床的利用において、上記コヒーレントレグレッション法に、さらに上記ＡＵＣ法を組み合わせることによって、従来の方法では計測が困難あるいは不可能であった、細い動脈の時間濃度曲線に関しても、パーシャルボリューム効果およびランダムなノイズの影響を排除し、しかも、未知の定数Ａを含まない測定値が得られる。
【０１４８】
なお、もちろんＡＵＣ法は、単独で従来の方法で計測された動脈に関する時間濃度曲線＜ｔ_ｋ，（ｖ´_ｋ（ｘ）−ｖ´_１（ｘ））＞に対しても適用でき、（ランダムなノイズやパーシャルボリューム効果の影響は排除できないものの、）未知であった定数Ａの値を決定した時間濃度曲線を構成できる。
【０１４９】
（上矢状静脈洞の時間濃度曲線を使った脳動脈の時間濃度曲線の補正（パーシャルボリューム効果の影響を抑圧））
パーシャルボリュームの影響を抑圧するために、このコヒーレントレグレッションに代えて又は併用して、上矢状静脈洞の時間濃度曲線Ｃｓｓｓ（ｔ）を使って脳動脈の時間濃度曲線Ｃａ（ｔ）を補正するようにしてもよい。
【０１５０】
まず、図５のステップＳ７において、図１２に示すように、ＣＴ画像上で上矢状静脈洞を取り囲むように大き目に上矢状静脈洞ＲＯＩを設定する。上矢状静脈洞は動脈に比べて大きく、また位置も比較的固定しているので、上矢状静脈洞ＲＯＩ設定は容易である。この大き目の上矢状静脈洞ＲＯＩには、複数の画素が含まれている。
【０１５１】
次に、上矢状静脈洞ＲＯＩの中の全ての画素がその全域に渡って上矢状静脈洞に含まれるように、上矢状静脈洞ＲＯＩを縮小処理する（Ｓ８）。縮小処理としては、例えば、まず、上矢状静脈洞ＲＯＩ内の画素各々について、しきい値処理（二値化）を実行し、ＲＯＩ内の二値マップ（“０”，“１”）を作成する。当該しきい値は、上矢状静脈洞の像を、その周辺の組織や骨の像と分離する値に設定される。“１”は上矢状静脈洞の像上の画素であることを表し、“０”は周辺の組織や骨の像上の画素であることを表している。この二値マップの各画素（中心画素）を、その近傍４または８画素の値に従って、置換する。中心画素が“１”であって、且つ、近傍４または８画素全てが“１”である場合のみ、中心画素の値を“１”に維持する。つまり、中心画素が“０”である場合はもちろん、たとえ“１”であったとしても、近傍４または８画素のなかの１つでも“０”を示しているときには、当該中心画素の値を“０”に置換する。従って、上矢状静脈洞ＲＯＩは、上矢状静脈洞の像の外形よりも、少なくとも１画素分縮小される。それにより縮小処理を受けた上矢状静脈洞ＲＯＩの中の全ての画素は、上矢状静脈洞像上の画素であるという条件を高い確度で実現され得る。
【０１５２】
また、この手法に代えて、時間濃度曲線の曲線下面積ＡＵＣを使って上矢状静脈洞ＲＯＩを修正するようにしてもよい。この場合、大き目のＲＯＩを探索範囲として、その中の画素各々について時間濃度曲線の曲線下面積ＡＵＣを計算する。造影効果により上矢状静脈洞像上の画素の曲線下面積ＡＵＣは、周辺画素のそれに比べて明らかに高値を示す。従って、この曲線下面積ＡＵＣに対してしきい値処理を実行することにより、ＲＯＩの中から上矢状静脈洞像上の画素だけを選別することができる。
【０１５３】
このようにしていずれかの手法又は両手法を併用してアンド条件によりピックアップされた上矢状静脈洞像上であるという確度の高い複数の画素に対して、各画素の時間濃度曲線が平均化され、上矢状静脈洞の時間濃度曲線Ｃｓｓｓ（ｔ）が作成される（Ｓ９）。
【０１５４】
ここで、ヨード造影剤は血液脳関門(Blood Brain Barrier)を通過しないので、原理的に、ヨード濃度は脳動脈と脳静脈とで変化しない、つまり、上矢状静脈洞の時間濃度曲線Ｃｓｓｓ（ｔ）の曲線下面積ＡＵＣは、Ｓ６で作成した脳動脈の時間濃度曲線Ｃａ（ｔ）の曲線下面積ＡＵＣに略等価になる。従って、図１３に示すように、上矢状静脈洞の時間濃度曲線Ｃｓｓｓ（ｔ）の曲線下面積ＡＵＣｓｓｓに対して、Ｓ６で作成した脳動脈の時間濃度曲線Ｃａ（ｔ）の曲線下面積ＡＵＣａが略等価になるように、Ｓ６で作成した脳動脈の時間濃度曲線Ｃａ（ｔ）の各時刻値に、ＡＵＣ（ｓｓｓ／ＡＵＣａ）を乗算することで補正する（Ｓ１０）。
【０１５５】
次に、以上のようにノイズ及びパーシャルボリューム効果が抑圧された図１４（ａ）に示す脳動脈の時間濃度曲線Ｃａ（ｔ）を使って、脳組織（毛細血管）の血流動態の様子を定量化する。そのためにまず、脳組織上の各画素について、図１４（ｂ）に示す時間濃度曲線Ｃｉ（ｔ）が作成される（Ｓ１１）。
【０１５６】
次に、Ｓ１２に示すように、左右エリアで別々の脳動脈時間濃度曲線Ｃａ（ｔ）を使って、画素ごとに、脳動脈時間濃度曲線Ｃａ（ｔ）を入力関数、脳組織の時間濃度曲線Ｃｉ（ｔ）を出力関数として、トレーサーが毛細血管を通過する過程の特徴を、伝達関数ＭＴＦとして求める。つまり、左エリアの脳組織の時間濃度曲線Ｃｉ（ｔ）に対しては同じ左エリアの脳動脈の時間濃度曲線Ｃａ（ｔ）を使い、また右エリアの脳組織の時間濃度曲線Ｃｉ（ｔ）に対しては同じ右エリアの脳動脈の時間濃度曲線Ｃａ（ｔ）を使って、伝達関数ＭＴＦを求める。また、上述したように脳動脈時間濃度曲線Ｃａ（ｔ）は、ＡＣＡ，ＭＣＡ，ＰＣＡごとに作成されるので、各Ｃａ（ｔ）ごとに伝達関数ＭＴＦの計算が繰り返される。
【０１５７】
ここでは、図１５に示すように、脳動脈の時間濃度曲線Ｃａ（ｔ）と、毛細血管の時間濃度曲線Ｃｉ（ｔ）との間で成り立つ理想的な関係を解析モデルとして用いてｂｏｘ−ＭＴＦ法を適用する。
【０１５８】
図１６にｂｏｘ−ＭＴＦ法の原理を示している。脳動脈の時間濃度曲線Ｃａ（ｔ）と、矩形関数で表される伝達関数ｂｏｘ−ＭＴＦとのコンボルーションＣｉ´（ｔ）と、Ｓ１１で作成した実測Ｃｉ（ｔ）との残差を評価し、その残差の二乗和を減少させるように伝達関数ｂｏｘ−ＭＴＦを修正する。この手続き(procedure)ルーチンを繰り返すことで、残差を最小化させる。
【０１５９】
残差を最小化する伝達関数ｂｏｘ−ＭＴＦに基づいて、図１５に示すように、ＣＢＰ、ＣＢＶ、ＭＴＴが計算され（Ｓ１３）、またＳ１２で最小化した残差の二乗和をＥｒｒとして出力する。厳密には、
ＣＢＰ＝ＣＢＰ
ＣＢＶ＝（１−Ｈｔ）／（１−ｂ＊Ｈｔ）＊ＣＢＶ'
ＭＴＴ＝（１−Ｈｔ）／（１−ｂ＊Ｈｔ）＊ＭＴＴ'
で補正される。ここで、Ｈｔは大血管のヘマトクリット値であり、ｂ＊Ｈｔは末梢血管のヘマトクリット値（一般的にはｂは０．７程度）である。
【０１６０】
なお、残差は、ｙ_ｉ(x)−ｆ(t_ｉ,x)で与えられる。ｙ_ｉ(x)は、時刻ｔ_ｉにおけるボクセルｘのスカラー値を表し、脳組織の時間濃度曲線に対応する。ｆ(t_ｉ,x)は、ボクセルｘのベクトル画素値に対してフィッティングされたモデルの時刻ｔ_ｉにおけるスカラー値を表し、伝達関数と脳動脈の時間濃度曲線とのコンボルーションに対応する。Ｅｒｒは、伝達関数を近似する際の残差の２乗和の平方根であり、例えば次式に示すように計算される。
【０１６１】
【数１４】

【０１６２】
ここで、Ｓは定数であり、例えばＳ＝Ｎ−ｐとする。ｐは、自由度、すなわち近似されるモデルｆに含まれるパラメータの数を表している。ｗ(t_ｉ)は時刻t_ｉにおける残差がＥｒｒに寄与する度合いを決める重み係数である。例えば、ｗ(t_ｉ)は、ｉに依存しないで、１等の固定値としても良いし、
ｗ(t_ｉ)＝αｅ^-ti ^２ ^/ ^β
として、｜ｔ_ｉ｜が大きいほど、重みｗが緩やかに小さくなるように構成し、時間が経つにつれて残差がＥｒｒにあまり寄与しなくなるようにしてもよい。
【０１６３】
こうしてｂｏｘ−ＭＴＦ法により得られた伝達関数から計算されたＣＢＰ、ＣＢＶ、ＭＴＴ、Ｅｒｒに対して、図１７に示すように、個々に出力範囲（適正範囲）が設定される（Ｓ１４）。それぞれ対応する出力範囲内の値を持つ画素については、その値が維持され、出力範囲外の値を持つ画素にはその値から、例えば表示上で黒に対応する値に置換される（Ｓ１５）。
【０１６４】
次に、図１８（ａ）乃至図１８（ｄ）に示すように、出力適正化を受けたＣＢＰ、ＣＢＶ、ＭＴＴ、Ｅｒｒからそれぞれのマップが生成される（Ｓ１６）。ＣＢＰ、ＣＢＶ、ＭＴＴ、Ｅｒｒのインデックスは、前大脳動脈ＡＣＡ，中大脳動脈ＭＣＡ，後大脳動脈ＰＣＡごと、しかも各スライスに対して個々に計算される。従って、マップは、１スライスでも、図１９に示すように、４×３＝１２枚になる。設定する脳動脈数を増加すれば、その増加数の４倍でマップは増える。このような多くのマップを総合的に評価することは現実的ではない。そこで、マップ枚数を減少させるために、マップを合成する（Ｓ１７）。
【０１６５】
図２０に示すように、合成方法としては、前大脳動脈ＡＣＡのＣＢＰマップと、中大脳動脈ＭＣＡのＣＢＰマップと、後大脳動脈ＰＣＡのＣＢＰマップとを、前大脳動脈ＡＣＡの残差Ｅｒｒ、中大脳動脈ＭＣＡの残差Ｅｒｒ、後大脳動脈ＰＣＡの残差Ｅｒｒに基づいて合成する。例えば前大脳動脈ＡＣＡの時間濃度曲線Ｃａ（ｔ）とその支配下にある脳組織の時間濃度曲線Ｃｉ（ｔ）とから伝達関数ＭＴＦを求めた場合、その残差Ｅｒｒは比較的少なく、逆に、支配下にない脳組織の時間濃度曲線Ｃｉ（ｔ）から伝達関数ＭＴＦを求めた場合、その残差Ｅｒｒは比較的多くなる。つまり、残差Ｅｒｒは各脳動脈の支配可能性を表している。
【０１６６】
従って、各画素ごとに、前大脳動脈ＡＣＡのＣＢＰ値と、中大脳動脈ＭＣＡのＣＢＰ値と、後大脳動脈ＰＣＡのＣＢＰ値との中から、最も残差Ｅｒｒの低い値に対応するＣＢＰ値をその画素の値として選択する。こうして合成されたマップは、前大脳動脈ＡＣＡ，中大脳動脈ＭＣＡ，後大脳動脈ＰＣＡの支配下にある可能性の高い脳組織のＣＢＰ値から構成される。他のインデックスＣＢＶ、ＭＴＴのマップ合成についても同様である。
【０１６７】
ここで、マップ合成について以下に詳細に説明する。動脈に対応する位置にある画素から得られる時間濃度曲線は動脈血中造影剤濃度を反映しており、これに上記したコヒーレントレグレッション法等を適用して正確な動脈血中造影剤濃度の時間濃度曲線を得ることができる。このような脳動脈の時間濃度曲線は動脈ごとに作成可能であり、それぞれ血行状態によって違いがある。特に脳血管障害等を起こしている症例に於いてはこの違いが著しい場合がある。例えばＫ箇所の動脈で得た脳動脈の時間濃度曲線をＡ_ｋ（t）（ｋ＝１,２，・・・，Ｋ）で表すことにする。
【０１６８】
ある組織の時間濃度曲線を、その組織を栄養している（支配している）動脈の脳動脈の時間濃度曲線と比較することによって、当該組織に於ける微小循環（毛細血管系の構造、機能）を反映するＣＢＰ等のインデックスを得ることが出来る。これらのインデックスは各部位（x,y,z）ごとに算出されるので、その値を画素値とする画像を構成することができ、このような画像が、インデックスマップである。例えばＲ種類（典型的には上述したようにＣＢＰ，ＣＢＶ，ＭＴＴ，Ｅｒｒの４種）のインデックスが得られる場合、Ｒ枚のマップが構成できる。このようにして作成されたＲ枚のマップは、各画素がベクトル値をとる１枚のマップ（ベクトル値マップ）と見なすことができる。すなわち、
Ｖ_ｋ（x,y,z）＝ <Ｐ_ｋ _, _１（x,y,z）, Ｐ_ｋ _, _２（x,y,z）, ... , Ｐ_ｋ _, _Ｒ（x,y,z）>
となる。
【０１６９】
例えばＣＢＰスタディでは、典型的には上述したようにＲ＝４とし、Ｐ_ｋ _, _１（x,y）はＣＢＰの値を、Ｐ_ｋ _, _２（x,y,z）はＣＢＶの値を、Ｐ_ｋ _, _３（x,y,z）はＭＴＴの値を、Ｐ_ｋ _, _４（x,y,z）は残差Ｅｒｒの値を表すように構成できる。
【０１７０】
部位（x,y,z）の内で、解析の対象となる臓器に対応していないことが明らかであるようなものは初めから解析の対象外とし、Ｐ_ｋ _,r（x,y）には解析対象外を示す特殊な値を代入するとよい（上記ステップＳ１４，Ｓ１５）。そのような値として、負で絶対値が大きい値を用いると便利である。或いは、ベクトルＶ_ｋ（x,y,z）に追加されるべきさらにもう一個の要素として、
Ｐ_ｋ _, _Ｒ ₊ _１（x,y,z）＝（（x,y,z）が解析対象外であれば０, さもなくば１）というマップを作っても良い。このようなマップは「マスク」と呼ばれる。
【０１７１】
このようなベクトル値マップＶ_ｋは、参照した脳動脈の時間濃度曲線Ａ_ｋごとに作成される。例えば、左右の中大脳動脈、前大脳動脈、後大脳動脈から脳動脈の時間濃度曲線を得たとするとＫ＝６、さらにおよび病変部周辺にある動脈数カ所から脳動脈の時間濃度曲線を得たとすると、Ｋ＝１０〜１５程度になる。
【０１７２】
このうち、右半球にある動脈から得られた脳動脈の時間濃度曲線は右半球に属する部位（x,y,z）の解析にだけ、また左半球にある動脈から得られた脳動脈の時間濃度曲線は左半球に属する部位（x,y,z）の解析にだけ、用いられるべきである。そこで、右半球と左半球の境目（正中線）を直線、曲線もしくは折れ線、或いは平面、曲面等として操作者が指定し、それぞれの半球ごとにマップを作るように構成するのが望ましい。しかしそれでもなお、片側の半球ごとにＫ＝３〜１０程度の数の脳動脈の時間濃度曲線が存在しうる。
【０１７３】
このように脳動脈の時間濃度曲線Ａ_ｋ（ｋ＝１,２，・・・，Ｋ）の数Ｋが大きい場合に、結果として得られるベクトル値マップＶ_ｋ（ｋ＝１,２，・・・，Ｋ）の枚数が多いために、観察するのに不便である。すなわち通常のグレースケール画像あるいはカラースケール画像として観察しようとすれば、一つのマップがＲ枚の画像から構成され、これがＫ組あるのだから、合計Ｋ×Ｒ枚の画像を比較しなくてはならない。さらに、どの部位がどの動脈によって栄養されているのかは必ずしも自明でなく、解剖学的知識を用いて、各部位ごとにどのマップＶ_ｋ（ｋ＝１,２，・・・，Ｋ）を観察するべきかを判断しなくてはならない。特に脳梗塞等の脳血管障害を生じている症例においては、組織を支配しているのがどの動脈かは、必ずしも解剖学的知識とは一致せず、異常な支配がしばしば見られる。これらの問題によって、ベクトル値マップの読影が難しいという問題点がある。
【０１７４】
この問題を解決するために、マップ合成を行う。つまり、残差マップを利用して、Ｋ個のベクトル値マップＶ_ｋ（ｋ＝１,２，・・・，Ｋ）を一つのベクトル値マップＶに集約する。例えば、Ｐ_ｋ _, _Ｒ（x,y,z）が残差マップである場合に、Ｖ（x,y,z）＝Ｖ_ｋ（x,y,z）ただしｋはｋ＝１,２，・・・，Ｋのうちで｜Ｐ_ｋ _, _Ｒ（x,y,z）｜が最小であるようなｋとする。
【０１７５】
また、各部位においてｋ＝１,２，・・・，Ｋのうちどれが採用されたかを示すためのマップ
Ｐ_０（x,y,z）＝（ｋ＝１,２，・・・，Ｋのうちで｜Ｐ_ｋ _, _Ｒ（x,y,z）｜が最小であるようなｋ）
を追加することもできる。
【０１７６】
この方式によれば、すなわち通常のグレースケール画像あるいはカラースケール画像として観察しようとするときＲ枚ないしＲ+１枚の画像を観察すればよい。
【０１７７】
この方式によれば、本来脳動脈の時間濃度曲線Ａ_ｊを使って算出されるべき部位（x,y,z）において誤ってＡ_ｋを使った算出結果が用いられる可能性がある。しかしながらこのような誤りが生じるには、Ｖ（x,y,z）の定義から明らかなように、｜Ｐ_ｋ _, _Ｒ（x,y,z）｜<｜Ｐ_j, _Ｒ（x,y,z）｜
となることが必要であり、このような関係は、Ａ_ｊとＡ_ｋが極めて類似している場合にしか生じない。このため、部位（x,y,z）においてはＶ_ｋ（x,y,z）とＶ_j（x,y,z）は元々類似していると考えられ、この誤りによって、Ｖ（x,y,z）の解釈に誤りが生じる可能性はほとんどない。
【０１７８】
実際にこの方法を適用すると、Ａ_ｊとＡ_ｋが極めて類似している場合にだけ、概ね一様であると思われる組織内において、部位ごとにＰ_０（x,y,z）＝ｋであったりＰ_０（x,y,z）＝jであったりすることが起こり、そのときＰ_ｋ _,r（x,y,z）≒Ｐ_j,r（x,y,z）（r＝１,２,...,Ｒ）であって、どちらを採用しても結果はほとんど違いがないことが観察される。
【０１７９】
逆に、特定の動脈に支配されている組織であって、それに対応する脳動脈の時間濃度曲線Ａ_ｋが他のカーブと似ていない場合には、本方式を用いることによって、当該組織中の部位（x,y,z）においてはほぼ確実に、しかも自動的にＶ_ｋ（x,y,z）が選択される。従って、上記Ｐ_０（x,y,z）を観察することによって、解剖学的知識なしに、どの組織がどの動脈の支配を受けているかを観察することができる。
【０１８０】
ここで、図６に戻る。図２１（ａ）乃至図２１（ｄ）に示すように、以上のように合成された、又は各脳動脈で単独のＣＢＰマップ、ＣＢＶマップ、ＭＴＴマップ、Ｅｒｒマップに対して、複数画素を含む関心領域ＲＯＩを設定し（Ｓ１８）、そのＲＯＩ内の画素値（ＣＢＰ値，ＣＢＶ値、ＭＴＴ値、Ｅｒｒ値）の平均値（ＣＢＰ平均値、ＣＢＶ平均値、ＭＴＴ平均値、Ｅｒｒ平均値）を計算し（Ｓ１９）、その平均値を診断材料とすることがある。この平均化に際して、上記ステップＳ１４でＣＢＰ、ＣＢＶ、ＭＴＴ、Ｅｒｒ各々に対して適正範囲を設定し、その範囲内の値を維持し、その範囲から外れた値は、例えば黒色表現に対応した最小値に置換したので、この置換した値を含めて平均化すると、その平均値には当然にして誤差が含まれてしまう。そのためこの平均化処理にあたっては、適正範囲内の値だけを選択して、または置換した値を除外して、平均化処理をすることが必要である。
【０１８１】
また、この平均化のための関心領域ＲＯＩの設定にあたっては、ＣＢＰマップ、ＣＢＶマップ、ＭＴＴマップ、Ｅｒｒマップのいずれかのマップ上で当該関心領域ＲＯＩを設定すれば、そのＲＯＩが他のマップにも共通で用いられるようになっており、、ＲＯＩ設定作業の簡素化を図り、また、同じＲＯＩに関する平均値（ＣＢＰ平均値、ＣＢＶ平均値、ＭＴＴ平均値、Ｅｒｒ平均値）の計算を可能としている。
【０１８２】
ここで、上述したように、ある脳動脈の時間濃度曲線に対するある組織の時間濃度曲線の最小残差Ｅｒｒは、その脳動脈がその組織をどの程度、支配しているか、つまりその脳動脈がその組織への血流供給をどの程度担っているかを表している。換言すると、その組織がその脳動脈からどの程度、従属しているか、つまりその脳組織がその脳動脈からどの程度血流供給を受けているかを表している。小さい残差Ｅｒｒに対応する脳動脈は、その画素の脳組織に対する支配可能性が高いことを表し、大きい残差Ｅｒｒに対応する脳動脈は、その画素の脳組織に対する支配可能性が低いことを表している。従って、残差Ｅｒｒから画素ごとに支配可能性の高い脳動脈をラベルで区別したマップ、つまり前大脳動脈ＡＣＡの支配可能性の高い領域と、中大脳動脈ＭＣＡの支配可能性の高い領域と、後大脳動脈ＰＣＡの支配可能性の高い領域とを区別した支配マップを生成することができる。
【０１８３】
図２２に示すように、脳の左右エリアそれぞれについて、前大脳動脈ＡＣＡの残差Ｅｒｒと、中大脳動脈ＭＣＡの残差Ｅｒｒと、後大脳動脈ＰＣＡの残差Ｅｒｒとを画素ごとに比較する。最も小さい残差を示す脳動脈（ＡＣＡ，ＭＣＡ又はＰＣＡ）が当該画素の脳組織を支配している可能性の最も高いことを表している。画素ごとに、最も支配可能性が高い、つまり残差Ｅｒｒが最も小さい値を示す脳動脈を特定する。特定した脳動脈に対応するラベルを各画素に与える。
【０１８４】
図２３に、生成した支配マップの例を示す。この支配マップは、ラベルを色や濃淡で区別して表示される。また、インデックスマップを任意のラベルでフィルタすることにより、図２４（ａ），図２４（ｂ），図２４（ｃ）に示すように、インデックスマップから、脳動脈（ＡＣＡ，ＭＣＡ又はＰＣＡ）ごとにその支配可能性の高い領域を抽出することができる。
【０１８５】
以上のように本実施形態によれば、コヒーレントフィルタ又はコヒーレントレグレッションにより、空間及び時間分解能の低下を抑えて、ノイズを抑制し、それによりＣＢＰスタディの解析精度を向上することができる。
【０１８６】
（変形例）
本発明は、上述した実施形態に限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で種々変形して実施することが可能である。さらに、上記実施形態には種々の段階が含まれており、開示される複数の構成要件における適宜な組み合わせにより種々の発明が抽出され得る。例えば、実施形態に示される全構成要件から幾つかの構成要件が削除されてもよい。
【０１８７】
【発明の効果】
本発明によれば、ＣＢＰスタディにより得られるマップの診断効率を向上することができる。
【図面の簡単な説明】
【図１】ＣＢＰスタディの原理説明図。
【図２】本発明の実施形態に係る脳組織内毛細血管の血流動態に関するインデックス演算装置の構成を示すブロック図。
【図３】本実施形態のコヒーレントフィルタによる画像処理の説明図。
【図４】本実施形態におけるコヒーレントフィルタによるノイズ抑制処理の流れを示すフローチャート。
【図５】本実施形態におけるインデックス演算処理全体の前半部のフローチャート。
【図６】本実施形態におけるインデックス演算処理全体の後半部のフローチャート。
【図７】図５のステップＳ３の分割線の一例を示す図。
【図８】図６のステップＳ１２の処理工数を減少させるために用いられるボロノイ図上の各脳動脈の勢力域を示す図。
【図９】図５のステップＳ４のＡＴマップ，ＰＴマップ，ＴＴマップを示す図。
【図１０】図５のステップＳ４のＡＴ，ＰＴ，ＴＴを示す図。
【図１１】図５のステップＳ６でスライス間で共通される脳動脈ＲＯＩを示す図。
【図１２】図５のステップＳ７で設定される上矢状静脈洞ＲＯＩを示す図。
【図１３】図５のステップＳ１０の脳動脈の時間濃度曲線の補正に関する補足図。
【図１４】図５のステップＳ１０、Ｓ１１で作成された脳動脈の時間濃度曲線Ｃａ（ｔ）と脳組織の時間濃度曲線Ｃｉ（ｔ）との一例を示す図。
【図１５】図６のステップＳ１２のｂｏｘ−ＭＴＦ法の原理説明図。
【図１６】図６のステップＳ１２のｂｏｘ−ＭＴＦ処理の説明図。
【図１７】図６のステップＳ１４の各インデックスの出力範囲設定画面の一例を示す図。
【図１８】図６のステップＳ１６で作成されたＣＢＰ，ＣＢＶ，ＭＴＴ，Ｅｒｒの各マップの一例を示す図。
【図１９】図６のステップＳ１６で脳動脈ごとに作成されたＣＢＰマップ，ＣＢＶマップ，ＭＴＴマップ，Ｅｒｒマップを一覧で表示した図。
【図２０】図６のステップＳ１７のマップ合成法を説明するための図。
【図２１】図６のステップＳ１９で計算された平均値の表示例を示す図。
【図２２】図６のステップＳ１７において、画素（局所組織）ごとに支配可能性の高い脳動脈（ＡＣＡ，ＭＣＡ，ＰＣＡ）をラベルにより区別したマップ（支配マップ）の生成方法を示す図。
【図２３】図２２の生成法により生成された支配マップの例を示す図。
【図２４】支配マップを使ってフィルタしたＣＢＰマップを示す図。
【符号の説明】
１０…ガントリ部、
２０…コンピュータ装置、
３０…画像処理装置、
１０１…Ｘ線管、
１０１ａ…高電圧発生装置、
１０２…Ｘ線検出器、
１０３…データ収集部、
１０４…前処理部、
１０５…メモリ部、
１０６…画像再構成部、
１０Ｍ…記憶装置、
１０７…画像表示部、
１０８…制御部、
１０９…入力部、
１１０…コヒーレントフィルタ処理部、
１１１…分散値推定部、
１１２…重み関数演算部、
１１３…画素値演算部（コヒーレントフィルタ部）、
１２０…ＣＢＰスタディ処理部、
１２１…ＲＯＩ設定支援部、
１２２…時間濃度曲線作成部、
１２３…脳動脈時間濃度曲線補正部、
１２４…ＭＴＦ処理部、
１２５…インデックス計算部、
１２６…マップ作成部、
１２７…マップ合成部。

Claims

造影剤を注入された被検体の特定部位に関する連続的な複数の画像から、前記特定部位内の動脈の設定を支援するために、前記画像を構成する複数の画素各々の時間濃度曲線の特性値を計算し、特性値マップを表示し、
前記設定された動脈に関する第１時間濃度曲線と、前記特定部位内の組織に関する第２時間濃度曲線を作成し、
前記動脈各々に対する前記組織の局所血流動態を表す伝達関数を、前記伝達関数と第１時間濃度曲線とのコンボルーションに対する前記第２時間濃度曲線の差が最小化するように計算し、
前記伝達関数から、前記動脈に対する前記局所血流動態に関するインデックスを計算して表示することを特徴とするインデックス計算方法。
前記特定部位は、脳であることを特徴とする請求項１のインデックス計算方法。
前記インデックスとして、脳組織の単位体積及び単位時間あたりの血流量、脳組織内の単位体積あたりの血液量、脳組織の血流平均通過時間の少なくとも１つが計算されることを特徴とする請求項２のインデックス計算方法。
前記伝達関数は、前記脳の左エリアと右エリアとで別々に計算されることを特徴とする請求項２記載のインデックス計算方法。
前記伝達関数は、前記脳の領域を左脳と右脳ごとに前大脳動脈（ＡＣＡ）、中大脳動脈（ＭＣＡ）、後大脳動脈（ＰＣＡ）各々上の母点としてボロノイ図法により分割されたＡＣＡ領域とＭＣＡ領域とＰＣＡ領域で別々に計算されることを特徴とする請求項２記載のインデックス計算方法。
前記画像から作成された前記静脈の時間濃度曲線の曲線下面積に、前記第１時間濃度曲線の時間濃度曲線の曲線下面積が略等価になるように、前記動脈に関する第１時間濃度曲線を補正することを特徴とする請求項２記載のインデックス計算方法。
前記静脈の時間濃度曲線は、前記静脈の領域内の複数画素に関する平均濃度から作成されることを特徴とする請求項６記載のインデックス計算方法。
前記特性値は、アピエランスタイム、ピークタイム又はトランジットタイムであることを特徴とする請求項１記載のインデックス計算方法。
前記インデックスの値が所定範囲から外れているとき、前記インデックスの値を黒、その近似色又は読影に使用しない色に対応する所定の値に置き換えることを特徴とする請求項１記載のインデックス計算方法。
造影剤を注入された被検体の特定部位に関する連続的な複数の画像から、前記特定部位内の動脈の設定を支援するために、前記画像を構成する複数の画素各々の時間濃度曲線の特性値を計算し、特性値マップを生成する特性値マップ生成部と、
前記設定された動脈に関する第１時間濃度曲線と、前記特定部位内の組織に関する第２時間濃度曲線を作成する時間濃度曲線作成部と、
前記動脈各々に対する前記組織の局所血流動態を表す伝達関数を、前記伝達関数と第１時間濃度曲線とのコンボルーションに対する前記第２時間濃度曲線の差が最小化するように計算する伝達関数計算部と、
前記伝達関数から、前記動脈に対する前記局所血流動態に関するインデックスを計算するインデックス計算部と、
前記インデックスのマップを作成するマップ作成部とを具備することを特徴とするインデックス計算装置。