(57)【要約】【課題】本発明の目的は、空間及び時間分解能の低下を抑えて、ノイズを抑制することにより、ＣＢＰスタディの解析精度を向上することにある。【解決手段】本発明は、脳血管透過性を持たない造影剤を注入された被検体の脳を撮影対象として連続的に取得した複数のＣＴ画像各々を構成する画素間の類似度を、複数のＣＴ画像にわたる各画素の時間濃度曲線に基づいて判定し、類似度に応じた重みにより画像ごとに画素を局所内で加重平均するコヒーレントフィルタ処理部１１０と、加重平均された複数のＣＴ画像から脳動脈画素の時間濃度曲線と脳組織画素の時間濃度曲線とを生成し、脳動脈画素の時間濃度曲線に対する脳組織画素の時間濃度曲線の伝達関数を計算し、伝達関数に基づいて脳組織内の毛細血管の血流動態に関するＣＢＰ，ＣＢＶ，ＭＴＴ，Ｅｒｒを演算するＣＢＰスタディ処理部１２０とを具備する。【選択図】図１
2008100121.tif 000002
【特許請求の範囲】
【請求項１】
脳血管透過性を持たない造影剤を注入された被検体の脳を撮影対象として連続的に取得した複数のＣＴ画像各々を構成する画素間の類似度を、前記複数のＣＴ画像にわたる各画素の時間濃度曲線に基づいて判定し、
前記類似度に応じた重みにより画像ごとに画素を局所内で加重平均し、
前記加重平均された複数のＣＴ画像から脳動脈画素の時間濃度曲線と脳組織画素の時間濃度曲線とを生成し、
前記脳動脈画素の時間濃度曲線に対する脳組織画素の時間濃度曲線の伝達関数を計算し、
前記伝達関数に基づいて、前記脳組織内の毛細血管の血流動態に関する複数種類のインデックスを演算することを特徴とする脳組織内毛細血管の血流動態に関するインデックス演算方法。
【請求項２】
時間軸に対する組織の濃度変化のインパルス応答が矩形関数で表されると仮定する解析モデルを用い、前記脳動脈画素の時間濃度曲線と前記伝達関数のコンボリューションにより得られる時間濃度曲線と、前記脳組織画素の時間濃度曲線との間の残差の指標が最小化するようにカーブフィッティングにより前記伝達関数を決定し、得られた伝達関数から前記インデックスを算出することを特徴とする請求項１記載のインデックス演算方法。
【請求項３】
前記インデックスとして、前記得られた伝達関数から、
ＣＢＰ：前記脳組織内の単位体積及び単位時間あたりの血流量  
ＣＢＶ：脳組織内の単位体積あたりの血液量  
ＭＴＴ：毛細血管の血液平均通過時間  
Ｅｒｒ：前記最小化された残差の二乗和  
とが計算されることを特徴とする請求項２記載のインデックス演算方法。
【請求項４】
前記類似度が比較的高い場合には比較的高い重み係数が与えられ、前記類似度が比較的低い場合には比較的低い重み係数が与えられることを特徴とする請求項１記載のインデックス演算方法。
【請求項５】
脳血管透過性を持たない造影剤を注入された被検体の脳を撮影対象として連続的に取得した複数のＣＴ画像各々を構成する画素間の類似度を、前記複数のＣＴ画像にわたる各画素の時間濃度曲線に基づいて判定する手段と、
前記類似度に応じた重みにより画像ごとに画素を局所内で加重平均する手段と、
前記加重平均された複数のＣＴ画像から脳動脈画素の時間濃度曲線と脳組織画素の時間濃度曲線とを生成する手段と、
前記脳動脈画素の時間濃度曲線に対する脳組織画素の時間濃度曲線の伝達関数をカーブフィッティングにより計算する手段と、
前記伝達関数に基づいて、前記脳組織内の毛細血管の血流動態に関する複数種類のインデックスを演算する手段とを具備することを特徴とする脳組織内毛細血管の血流動態に関するインデックス演算装置。
【請求項６】
脳組織内毛細血管の血流動態に関するインデックスをコンピュータに演算させるためのプログラムコードを記録したコンピュータ読み取り可能な記憶媒体であって、前記プログラムコードは、
脳血管透過性を持たない造影剤を注入された被検体の脳を撮影対象として連続的に取得した複数のＣＴ画像各々を構成する画素間の類似度を、前記複数のＣＴ画像にわたる各画素の時間濃度曲線に基づいて判定する手段と、
前記類似度に応じた重みにより画像ごとに画素を局所内で加重平均する手段と、
前記加重平均された複数のＣＴ画像から脳動脈画素の時間濃度曲線と脳組織画素の時間濃度曲線とを生成する手段と、
前記脳動脈画素の時間濃度曲線に対する脳組織画素の時間濃度曲線の伝達関数をカーブフィッティングにより計算する手段と、
前記伝達関数に基づいて、前記脳組織内の毛細血管の血流動態に関する複数種類のインデックスを演算する手段とを具備することを特徴とする記憶媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、脳組織内毛細血管の血流動態に関するインデックスの演算方法及び演算装置に関する。
【背景技術】
【０００２】
従来からＸ線ＣＴ検査では単純ＣＴ像から形態情報を、造影ＣＴによるダイナミックスキャンで病巣の周りの血流の動態情報をそれぞれ視覚情報として得ることができる。近年、マルチスライスＣＴによる高速スキャンが可能になり益々造影ＣＴのダイナミックスキャンの活用範囲が拡大していくものと考えられる。
【０００３】
その一つの方向性として、脳組織内毛細血管の血流動態に関するインデックスを演算するためのＣＢＰスタディと呼ばれる方法がある。ＣＢＰスタディでは、”毛細血管を通過する血流”の動態を定量的に表すＣＢＰ、ＣＢＶ、ＭＴＴ、Ｅｒｒ等のインデックスを求め、またこれらインデックスのマップを出力する。
【０００４】
ＣＢＰは、脳組織内の単位体積及び単位時間あたりの血流量［ml/100ml/miｎ］を表し、ＣＢＶは、脳組織内の単位体積あたりの血液量［ml/100ml］、ＭＴＴは毛細血管の血液平均通過時間［秒］を表している。
【０００５】
これら脳組織中の毛細血管の血流動態を定量的に表しているインデックスＣＢＰ、ＣＢＶ、ＭＴＴは、脳虚血が発症してからの経過時間情報ととともに、虚血性脳血管障害の病体鑑別、毛細血管の拡大の有無、血流速などの評価のための有益な情報として期待されている。例えば、一般に虚血性の脳血管障害では、提供する動脈の血圧が低下し、その血管内の血流速の低下が見られる。その結果、ＣＢＶは一定でも、ＭＴＴが延長し、ＣＢＰは低下する。また、脳梗塞超急性期では、血圧低下による血流速の低下を補うために、毛細血管を拡張させ、血流速を増加させることにより、血流量ＣＢＰの低下を抑制しようとする働き（オートレギュレーション）がある。従って、ＭＴＴが延長することにより、ＣＢＰが低下しても、ＣＢＶが増加していれば、毛細血管の再開通の可能性を示唆する情報となる。
【０００６】
ＣＢＰスタディではトレーサーとして脳血管透過性を持たない造影剤、例えばヨード造影剤が使用される。ヨード造影剤は例えばインジェクターにより肘静脈から注入される。インジェクターにより静注されたヨード造影剤は、心臓、肺を経由して、脳動脈から流れ込む。そして、脳動脈から、脳組織内の毛細血管を経て、脳静脈へと流れ出ていく。このとき、ヨード造影剤は正常な脳組織内の毛細血管では造影剤は血管外へ漏れ出ることなく通過する。図１はこの様子を模式的に示している。
【０００７】
造影剤の通過の様子をダイナミックＣＴで撮影して、その連続画像から、脳動脈の画素の時間濃度曲線Ｃａ（ｔ）、毛細血管を含む脳組織の画素の時間濃度曲線Ｃｉ（ｔ）、脳静脈の画素の時間濃度曲線Ｃｓｓｓ（ｔ）をそれぞれ測定する。
【０００８】
ここで、ＣＢＰスタディでは、造影剤の血中濃度について脳組織に近い脳血管の血中濃度Ｃａ（ｔ）と毛細血管の血中濃度Ｃｉ（ｔ）との間で成り立つ理想的な関係を解析モデルとしており、つまり脳組織に入る直前の血管から造影剤を注入した場合、毛細血管を含む脳組織単位体積（１画素）内の時間濃度曲線は立ち上がりが垂直で、しばらくは一定の値を維持し、その後、急勾配で立ち下がる形になる。これを矩形関数で近似する（ｂｏｘ−ＭＴＦ法：ｂｏｘ−Ｍｏｄｕｌａｔｉｏｎ  Ｔｒａｎｓｆｅｒ  Ｆｕｎｃｔｉｏｎ  ｍｅｔｈｏｄ）。
【０００９】
脳動脈血中時間濃度曲線Ｃａ（ｔ）を入力関数、脳組織の時間濃度曲線Ｃｉ（ｔ）を出力関数として、トレーサーが毛細血管を通過する過程の特徴を、矩形の伝達関数として求めることができる。
【００１０】
このようなＣＢＰスタディにおける第１の問題は次の通りである。肘静脈に、ボーラスインジェクションを行った場合、ＣＴで観察される造影効果は、血液のＣＴ値（造影されない時数十ＨＵ）が最大数百ＨＵに上昇する。しかし、脳血流を有効に解析するためには造影効果を高々数パーセント以内の誤差で計測できなくてはならない。すなわち、血液の造影効果（ＣＴ値の上昇）が２０〜４０ＨＵ程度であってもこれを検出できる必要がある。
【００１１】
単位体積の脳組織中に占める毛細血管の体積比率は高々３〜４パーセント程度である。従って、血液のＣＴ値が２０〜４０ＨＵ上昇した場合、脳組織の平均ＣＴ値は、０．５〜１．５ＨＵ程度上昇するに過ぎない。
【００１２】
ＣＴ画像ではノイズの標準偏差（ｓｄ）は、照射Ｘ線量の平方根に反比例し、典型的な照射条件においてｓｄは例えば５〜１０ＨＵ程度である。従って、０．５ＨＵの造影効果を検出するためには、Ｘ線量を１０〜１００倍程度増やさねばならず、これは患者の被曝線量が著しく大きくなることを意味する。さらに、ダイナミックＣＴにおいては同一箇所を数十回に渡って撮影するのであるから、撮影箇所に於ける皮膚の被曝は通常の数百〜数千倍に至ることになり、炎症・脱毛・壊死・発癌等の放射線障害が生じる恐れがある。
【００１３】
むしろダイナミックＣＴにおいてはＸ線量を通常の撮影よりも減らさなくてはならない。一般に、１スキャン当たりのＸ線量を例えば通常の１／２〜１／１０程度に減じることが行われる。これによって、通常の１回のＣＴ撮影に比べて数倍〜２０倍程度のＸ線被曝に留めることができ、これは放射線障害を生じない程度である。しかし、このようなＸ線量を低減したＣＴ画像において、ｓｄは例えば１５〜２０ＨＵ程度であり、０．５〜１．５ＨＵ程度の造影効果は到底検出できない。
【００１４】
そこで、画像のノイズ成分を抑制することが、ＣＢＰスタディでは最も重要な課題の１つである。そのため、１）スライス厚を厚くする、２）画素束ね、３）画像の平滑化、が一般的に取りうる方策である。しかしこれらには以下のような問題点がある。
【００１５】
“スライス厚を厚くする”ために、撮影時にスライス厚を厚く設定するか、連続する薄いスライスの画像数枚を平均して厚いスライスの画像を生成する。スライス厚に比例して画素当たりのＸ線量が増えるため、画像ノイズのｓｄは、スライス厚の平方根に反比例して小さくなる。しかしながら、スライス厚を厚くすることによって、パーシャルボリューム効果が生じ、すなわち１個の画素が、一様な脳組織を表して居らず、複数の組織（白質、灰白質、血管、脳溝、脳室など）の平均的なＣＴ値を表すことになる確率が大きくなり、解析結果として得られる脳血流量等の値の誤差が大きくなる。
【００１６】
特に血管の影響を含む画素は、正常な解析が不可能である。このためスライス厚を厚くすると、不正確で、しかも解析不可能な画素を沢山含む非常に品質の悪い結果しか得られなくなる。
【００１７】
“画素束ね”では、空間解像度が或る程度犠牲になる。例えば一辺がｎ個の画素からなる正方形の領域（ｎ×ｎ個の画素を含む）の平均値を求め、これをその正方形全体の平均ＣＴ値とし、このような正方形を画素とみなし、これを敷き詰めて「画素束ね画像」を構成する。もとの画像が例えば一辺５１２個の画素からなる（５１２×５１２個の画素を含む）とし、ｎ＝２とすれば、「画素束ね画像」は一辺（５１２／２）個の画素から構成される（２５６×２５６個の画素を含む）画像となる。この方法によれば、ノイズは、ｎに反比例して減少させることが可能である。さらに、解析対象となる画素の数が１／（ｎ×ｎ）倍になるため、計算量も小さくなるという利点がある。
【００１８】
しかしながら、ｎを大きくすると、空間解像度が低下し、それに伴ってパーシャルボリューム効果が生じ、すなわち１個の画素が、一様な脳組織を表して居らず、複数の組織（白質、灰白質、血管、脳溝、脳室など）の平均的なＣＴ値を表すことになる確率が大きくなり、解析結果として得られる脳血流量等の値の誤差が大きくなる。特に血管の影響を含む画素は、正常な解析が不可能である。このため、ｎを大きくすると、空間解像度が低く、不正確で、しかも解析不可能な画素を沢山含む非常に品質の悪い結果しか得られなくなる。このため、実用上は、ｎ＝２〜４程度が限界であり、これだけでは十分なノイズ抑制効果が得られない。
【００１９】
また、画像の平滑化、すなわち１枚のＣＴ画像ごとに、２次元の空間フィルタを作用させて平滑化を行う方法を用いると、十分なノイズ抑制効果と引き換えに、空間解像度が著しく損なわれる。特に、太い血管（動脈・静脈）が存在する箇所に近接している画素には、太い血管において生じた造影効果の影響が及ぶことになり、これらの画素の時間濃度曲線は正しくなくなってしまう。従ってごく軽度の平滑化を行うに留めねばならない。ここで、ごく軽度の平滑化を行うに際して重要なのは、画像フィルタのサイズをごく小さくする事、例えば、３×３程度に設定することである。３×３の平滑化フィルタを用いて最大の画像ノイズ抑制効果を得ようとすると、その上限は、ノイズｓｄを１／３に低減することであり、それ以上にノイズを抑制するのは不可能である。従って十分なノイズ抑制効果は得られない。
【００２０】
一方、時間的平滑化、すなわち各画素について得られた時間濃度曲線を曲線とみなして、これを１次元フィルターで平滑化する手法を用いると、十分なノイズ抑制効果を得ようとすると時間分解能を著しく損なう。元来、ＣＢＰスタディでダイナミックＣＴを行うのは短いサンプリング周期で撮影を行うことによって高い時間分解能を得て、時間濃度曲線の僅かな変化（特に生理学的構造に起因する平滑化効果がどの程度生じているか）を精密に計測することが目的であるから、時間的平滑化は全く適当でない。
【発明の開示】
【発明が解決しようとする課題】
【００２１】
本発明の目的は、空間及び時間分解能の低下を抑えて、ノイズを抑制することにより、ＣＢＰスタディの解析精度を向上することにある。
【課題を解決するための手段】
【００２２】
本発明は、脳血管透過性を持たない造影剤を注入された被検体の脳を撮影対象として連続的に取得した複数のＣＴ画像各々を構成する画素間の類似度を、前記複数のＣＴ画像にわたる各画素の時間濃度曲線に基づいて判定し、前記類似度に応じた重みにより画像ごとに画素を局所内で加重平均し、前記加重平均された複数のＣＴ画像から脳動脈画素の時間濃度曲線と脳組織画素の時間濃度曲線とを生成し、前記脳動脈画素の時間濃度曲線に対する脳組織画素の時間濃度曲線の伝達関数を計算し、前記伝達関数に基づいて、前記脳組織内の毛細血管の血流動態に関する複数種類のインデックスを演算する。
【発明の効果】
【００２３】
本発明によれば、コヒーレントフィルタ又はコヒーレントレグレッションを併用することにより、空間及び時間分解能の低下を抑えて、ノイズを抑制し、それによりＣＢＰスタディの解析精度を向上することができる。
【発明を実施するための最良の形態】
【００２４】
以下、図面を参照して本発明を好ましい実施形態により説明する。  
本実施形態の特徴としては、コヒーレントフィルタを用いて空間及び時間分解能の低下を抑えながら、ノイズを効果的に抑制することにより、脳組織の毛細血管の血流動態を定量的に表すＣＢＰ等のインデックスを高精度に計算することにある。
【００２５】
（装置構成）  
図２には、本実施形態に係るＸ線ＣＴ装置の構成を示している。Ｘ線ＣＴ装置は、ガントリ部１０とコンピュータ装置２０とから構成される。ガントリ部１０は、Ｘ線管１０１、高電圧発生装置１０１ａ、Ｘ線検出器１０２、データ収集部１０３（ＤＡＳ；Ｄａｔａ  Ａｑｕｉｓｉｔｉｏｎ  Ｓｙｓｔｅｍ）とを有する。Ｘ線管１０１とＸ線検出器１０２とは、高速で且つ連続的に回転する図示しない回転リングに被検体Ｐを挟んで互いに対向する位置に搭載される。
【００２６】
コンピュータ装置２０は、画像処理装置３０と、画像表示部１０７と、入力部１０９とから構成される。画像処理装置３０は、制御部１０８を中枢として、データ収集部１０３から出力される生データを補正処理等を経て投影データに変換する前処理部１０４、投影データを記憶するメモリ部１０５、投影データからＣＴ画像データを再構成する画像再構成部１０６、ＣＴ画像データを保管する記憶装置１０Ｍ、ＣＴ画像データに対してコヒーレントフィルタ処理を実行するコヒーレントフィルタ処理部１１０、及びコヒーレントフィルタ処理を受けたＣＴ画像データを使ってＣＢＰスタディ処理を実行するＣＢＰスタディ処理部１２０とから構成される。
【００２７】
コヒーレントフィルタ処理部１１０は、分散値推定部１１１、重み関数演算部１１２、画素値演算部（コヒーレントフィルタ部）１１３とから構成される。これら分散値推定部１１１、重み関数演算部１１２、画素値演算部１１３の機能については後述するコヒーレントフィルタ処理の詳細説明の中で説明する。
【００２８】
ＣＢＰスタディ処理部１２０は、ＲＯＩ設定支援部１２１、時間濃度曲線作成部１２２、脳動脈時間濃度曲線補正部１２３、ＭＴＦ処理部１２４、インデックス計算部１２５、マップ作成部１２６、マップ合成部１２７から構成される。
【００２９】
ＲＯＩ設定支援部１２１は、ＣＴ画像上に脳動脈や脳静脈に対して関心領域ＲＯＩを設定する作業を支援するための情報（脳動脈ＲＯＩのためのＡＴマップ、ＰＴマップ、ＴＴマップ等）を作成し提供する。
【００３０】
なお、脳動脈ＲＯＩは、前大脳動脈（ＡＣＡ）、中大脳動脈（ＭＣＡ）、後大脳動脈（ＰＣＡ）を対象として、左脳、右脳それぞれの領域に個別に設定される。従って、左右に３個ずつ、合計６個の脳動脈ＲＯＩが設定される。また、脳動脈の時間濃度曲線Ｃａ（ｔ）を補正するために、他の時間濃度曲線Ｃｓｓｓ（ｔ）が利用される。この時間濃度曲線Ｃｓｓｓ（ｔ）は、パーシャルボリュームを含まない画素が存在するのに充分に太い血管上に設定されたＲＯＩに関して作成される。例えば、脳血管の中で最も太い上矢状静脈洞に設定することが好ましい。
【００３１】
時間濃度曲線作成部１２２は、記憶装置１０Ｍに記憶されているダイナミックＣＴ画像データ（時間的に連続した複数枚の画像データ）から脳動脈、脳静脈及び脳組織（毛細血管）に関する時間濃度曲線を作成する。なお、脳動脈の時間濃度曲線Ｃａ（ｔ）は、設定された６つの脳動脈ＲＯＩに関して個々に作成される。脳静脈の時間濃度曲線Ｃｓｓｓ（ｔ）は、上矢状静脈洞に設定された脳静脈ＲＯＩに関して作成される。また、脳組織の時間濃度曲線Ｃｉ（ｔ）は、脳組織上の全画素を対象として画素ごとに作成される。
【００３２】
脳動脈時間濃度曲線補正部１２３は、ノイズやパーシャルボリューム効果の影響を除去するために、脳動脈の時間濃度曲線Ｃａ（ｔ）を、上矢状静脈洞の時間濃度曲線Ｃｓｓｓ（ｔ）に基づいて補正する。この補正方法については後述する。ＭＴＦ処理部１２４は、補正された脳動脈の時間濃度曲線Ｃａ（ｔ）と、脳組織の時間濃度曲線Ｃｉ（ｔ）とに基づいて、ｂｏｘ−ＭＴＦ法により、伝達関数ＭＴＦを、脳組織上の全画素を対象として画素ごとに計算する。
【００３３】
インデックス計算部１２５は、計算された伝達関数ＭＴＦから脳組織の血流動態を表すインデックス（ＣＢＰ、ＣＢＶ、ＭＴＴ、Ｅｒｒ）を、脳組織上の全画素を対象として画素ごと計算する。マップ作成部１２６は、計算されたインデックス各々のマップを、脳動脈（ＡＣＡ，ＭＣＡ，ＰＣＡ）ごとに生成する。このマップは、各スライスに関して、インデックスの種類（＝４）×脳動脈の種類の数（ＡＣＡ，ＭＣＡ，ＰＣＡの３つ）＝１２種類作成される。マルチスライスでは、そのスライス数倍の種類のマップが作成される。この膨大な枚数のマップの枚数を合成処理により減らして診断効率を向上させるためにマップ合成部１２７が設けられている。
【００３４】
以下に、コヒーレントフィルタ処理とＣＢＰスタディ処理について順番に説明する。
【００３５】
コヒーレントフィルタの原理について簡単に説明すると、近傍の３×３等の局所内画素を加重平均し、その加重平均値を局所中心画素の値とすることを基本として、周辺画素各々の重みを中心画素と周辺画素との間の類似度に従って変えることを特徴としたものである。ここで言う類似度とは、画素間で、解剖学的に近い組織、具体的には同じ脳動脈の支配下にある毛細血管どうしである可能性の度合いを示す指標であり、この類似度が高い画素に対しては高い重みを与え、逆に類似度が低い画素に対してはゼロに近い低い重みを与えることにより、ノイズ抑制を果たしながらも、空間分解能の低下を抑制することを可能としている。ここで類似度の計算が重要になるが、本実施形態では、脳血管透過性を持たない造影剤、例えばヨード造影剤を注入（静注）された被検体の脳を撮影対象として連続的に取得した複数のＣＴ画像（ダイナミックＣＴ画像）を用いて、各画素の時間濃度曲線の比較により類似度を計算する。そのため類似度の確からしさは、サンプリング周波数、つまり単位時間あたりの画像枚数と、サンプリング数、つまり全画像枚数とに依存して決まる。そこでスキャン間隔を例えば０．５秒に短縮することが効果的である。
【００３６】
（コヒーレントフィルタ）  
（コヒーレントフィルタの一般的説明）
（画素値ｖ（ｘ））  
一般に、カメラ等の撮像手段を介して取得されたデジタル画像は、複数の画素（ｐｉｘｅｌ）から構成されている（あるいは、当該画像をそのような画素の集合として考えることができる。）。以下の説明では、当該画素の位置をベクトルｘ（すなわち座標値のベクトル）として表し、画素ｘが有する値（例えば濃淡を表わす数値、ＣＴ値ＨＵ）をＫ次元ベクトルとして表す。２次元画像の場合、画素ｘとは画像上における位置を表す座標値（ｘ、ｙ）を示す２次元ベクトルである。ある画素ｘについて定義される「画素値ｖ（ｘ）」を、  
ｖ（ｘ）＝（ｖ１（ｘ），ｖ２（ｘ），…，ｖＫ（ｘ））  …  （１）
と表記する。なお、この（１）式の右辺における、ｖ１（ｘ），ｖ２（ｘ），…，ｖＫ（ｘ）それぞれを、以下では、画素ｘについての「スカラー値」と呼ぶことにする。
【００３７】
例えば、画像が「カラー画像」であるとき、各画素が、それぞれ三原色（赤，緑，青）の明るさ（スカラー値）を有することから、これら各画素の画素値ｖ（ｘ）は、その次元がＫ＝３のベクトルであると考えることができる（上記（１）式の右辺各項で、その添え字が例えば「赤」，「緑」及び「青」である場合を想定されたい。また例えば、画像がＫ枚の静止画像から構成される動画像であって、第ｎ番目の画像の各画素はスカラー値ｖｎ（ｘ）を持つという場合には、Ｋ枚の静止画像上、共通の同一点（同一座標）の画素ｘの持つ画素値（スカラー値）を並べて構成される、Ｋ次元ベクトル値ｖｎ（ｘ）＝（ｖ１（ｘ），ｖ２（ｘ），…，ｖＫ（ｘ））が以下で述べるベクトル値としての画素値である。
【００３８】
（類似度（適合度ないし危険率）ｐ（ｘ，ｙ）と重みｗ（ｐ（ｘ，ｙ）））  
上記画素ｘに対して、適当な画素の集合Ｎ（ｘ）を考える（この集合Ｎ（ｘ）は画素ｘを含んでよい。）。次に、Ｎ（ｘ）の要素であるそれぞれの画素ｙと、前記画素ｘとの間で、重みｗ（ｐ（ｘ，ｙ））を考える。この重みｗ（ｐ（ｘ，ｙ））は、次に記す性質を有する。
【００３９】
（類似度ｐ（ｘ，ｙ））  
まず、重みｗ（ｐ（ｘ，ｙ））の値を左右する関数ｐ（ｘ，ｙ）の意味について述べる。このｐ（ｘ，ｙ）は、本実施形態にいう「類似度」を定量化する手段であり、一般的にいえば、画素ｘと画素ｙ∈Ｎ（ｘ）とが、何らかの意味でどの程度類似しているか（例えば、両画素ｘ及びｙの上記画素値ｖ（ｘ）及びｖ（ｙ）間に認められる統計的差異の程度）、を示す具体的数値を与える。
【００４０】
より具体的には例えば、ｐ（ｘ，ｙ）が小さな値を与えるときには、画素ｘと画素ｙとが、その画素値ｖ（ｘ）及びｖ（ｙ）間に「統計的に有意な差がなく（つまり類似度が高い）」、類似である可能性が高いと判断され、ｐ（ｘ，ｙ）が大きな値を与えるときには、「統計的に有意な差があり（つまり類似度が低い）」、の如く判断されるということである。
【００４１】
ところで、画素値ｖ（ｘ）及びｖ（ｙ）（ないしスカラー値ｖ１（ｘ），…，ｖＫ（ｘ）及びｖ１（ｙ），…，ｖＫ（ｙ））には、必ずノイズが含まれていると考えなければならない。例えば、画像がＣＣＤ撮像素子により取得された場合を考えると、それを構成する各画素については、素子内の暗電流や外界から入射する光量の不規則変動に起因するノイズ等が存在する。
【００４２】
このようなノイズは、一般に、全画素についてまちまちな値をとるため、画素ｘと画素ｙとが、仮に（外界における）同一物体を反映したものである場合であっても、実際に観測される画像上では、同一の値を持たないことがある。このことを逆にいえば、いずれも同一物体を反映した画素ｘと画素ｙにおいて、それぞれのノイズを除去した状況を仮に想定すれば、これらは該同一物体を表象するものとして画像上に表示され（＝そのように認識され）るし、また、両者は本来同一の（あるいはごく近い）画素値を有する。
【００４３】
そこで、上述したノイズの性質を踏まえ、上記のｐ（ｘ，ｙ）に関し、統計的検定法でよく知られている「帰無仮説」の概念を用いると、このｐ（ｘ，ｙ）については、具体的に次のように言うことができる。すなわち、帰無仮説Ｈ「画素ｘと画素ｙとはそれぞれのノイズを除去した場合に同一の画素値を有する」言いかえれば「ｖ（ｘ）＝ｖ（ｙ）、ただし、両画素のノイズに起因する差異を除く」を立てる（つまり、このような命題が成立する場合、「両画素ｘ及びｙとの間の類似度が高い（適合度が大きい）」と考える。）と、関数ｐ（ｘ，ｙ）は、この仮説Ｈを棄却する場合の危険率（あるいは、有意水準）であるということができる（この場合、ｐ（ｘ，ｙ）は、その値域が[０，１]であるような関数として定義される（ｐ（ｘ，ｙ）∈[０，１]）。）。
【００４４】
したがって、危険率ｐ（ｘ，ｙ）が大きい場合、すなわち棄却が誤りである危険性が大きい場合には上記仮説Ｈを満たす可能性が高いといえ、逆に小さい場合、すなわち棄却が誤りである危険性が小さい場合には仮説Ｈを満たさない可能性が高いということができる（なお、統計的検定における周知事項ではあるが、仮説Ｈが「棄却」されないといっても、それが「真」であることを意味するわけではない。この場合、仮説Ｈが示す命題が、否定し得ないことを意味するに過ぎない。）。
【００４５】
（重みｗ（ｐ（ｘ，ｙ）））  
さて、重みｗ（ｐ（ｘ，ｙ））は、その表され方から明らかな通り、上記したような危険率ｐ（ｘ，ｙ）の関数（より一般には、適合度の関数（適合度をρ（ｘ，ｙ）とすれば、ｗ（ρ（ｘ，ｙ））となるように構成できる）であり、また、この重みｗ（ｐ（ｘ，ｙ））を求めるため、ｘ及びｙの組み合わせそれぞれについて求められた危険率ｐ（ｘ，ｙ）に作用させる重み関数ｗは、一般的にいうと、上記「棄却」を具現化する作用を有するものである。具体的には、危険率ｐ（ｘ，ｙ）が大きい場合には重み関数ｗの値、すなわち重みｗ（ｐ（ｘ，ｙ））が大きな正の値をとり、その逆の場合には小さな正の値（又は“０”）をとる、等というように調整されている（重み関数ｗの具体的形式については後述する。）。つまり、重みｗ（ｐ（ｘ，ｙ））は、画素ｘと画素ｙとが、上記仮説Ｈに示される命題を満たすらしい場合には、大きい値をとり、その逆の場合には小さい値をとる。一例として特に、ｗのとりうる値が”０”かまたは”０”でない一定値の２通りしかないように構成してもよい。
【００４６】
なお、以上までに述べた仮説Ｈ、危険率ｐ（ｘ，ｙ）、重みｗ（ｐ（ｘ，ｙ））間の関係をまとめると、帰無仮説Ｈが正しい可能性が高いとき、類似度ｐも高くなり、その画素に与える重みｗを高くし、一方、帰無仮説Ｈが正しい可能性が低いとき、類似度ｐも低くなり、その画素に与える重みｗを低くする。このように加重平均値への寄与度（重み）を類似度に応じて変えることにより、分解能の低下を抑えながら、ノイズを効果的に抑制することが可能となる。また、重み関数ｗ（ｔ）は、より一般に、「ｔ∈[０，１］で定義される非負の単調増加関数」ということができ、また、該ｗ（ｔ）の満たすべき性質は、少なくともそのようであればよい。
【００４７】
（コヒーレントフィルタ処理）  
以上までの説明により、「コヒーレントフィルタ」は次のように導かれる。すなわちまず、画像を構成するある画素ｘに対し、集合Ｎ（ｘ）の要素たる画素ｙのすべてについて上記した重みｗ（ｐ（ｘ，ｙ））を計算する。次に、これら複数の重みｗ（ｐ（ｘ，ｙ））を用いて、当該画素ｘを構成する新たなスカラー値ｖ´ｋ（ｘ）を、以下の（２）式で計算する。すなわち、
【数１】
2008100121.tif 000003
【００４８】
ただし、ｋ＝１，２，…，Ｋである。そして、この式で求められたｖ´ｋ（ｘ）を用いて、当該画素ｘの変換後の画素値（新たな画素値）ｖ´（ｘ）を、  
ｖ´（ｘ）＝（ｖ´１（ｘ），ｖ´２（ｘ），…，ｖ´Ｋ（ｘ））  …  （３）  
として構成する。
【００４９】
ここに、上記（１）式で表される、画素値ｖ（ｙ）＝（ｖ１（ｙ），ｖ２（ｙ），…，ｖＫ（ｙ））（ｙ＝ｘである場合を含む。）を、ｖ´（ｘ）＝（ｖ´１（ｘ），ｖ´２（ｘ），…，ｖ´Ｋ（ｘ））に変換するフィルタが、「コヒーレントフィルタ」の形式である。これはその表式から明らかな通り、画素値を構成するスカラー値ｖｋ（ｙ）の重み付け平均値を表している。
【００５０】
このような処理は、以下のような結果をもたらす。すなわち、画素値ｖ´（ｘ）は、画素ｘとノイズを除いて同一の画素値をとることが確からしい（＝上記仮説Ｈの命題を満たす可能性が高い）画素ｙを重視した重み付け平均値ｖ´ｋ（ｘ）から構成されたベクトルを表すこととなる。また、このような画素ｙが十分な数存在するならば、画素値ｖ´（ｘ）は、画素ｘが本来有すべきその真値から外れることなく、上記したような平均化の作用によりノイズのみを抑制した値を有することとなる。
【００５１】
なお、危険率ｐ（ｘ，ｙ）が小さく、したがって、帰無仮説Ｈが「棄却」され、重みｗ（ｐ（ｘ，ｙ））が小さくなるような場合であっても、上記記述からもわかる通り、必ずしもこれを完全に「棄却」するとは限らない。このようなことは、後述する重み関数ｗの具体的形式に依存するところであるが、危険率ｐ（ｘ，ｙ）が“０”（＝０％）に近いような場合でも、ｗ（ｐ（ｘ，ｙ））≠０（ただし、ｐ（ｘ，ｙ）が“１”に近い場合に比べて、より小さな正の値ではある。）としてよい（なお、ｐ（ｘ，ｙ）＝１である場合とは、後述するように、ｖ（ｘ）＝ｖ（ｙ）のときである。）。
【００５２】
すなわち、完全な棄却ということではなく、小さな寄与は認めてよいということである（なおこのような場合に、ｗ（ｐ（ｘ，ｙ））＝０とするのであれば、完全な棄却を行うのと同義である。
【００５３】
このような処理は、一般的に次のように言える。すなわち、ある画像を構成する（複数の）画素ｘが存在するとき、この画素ｘとある任意の画素ｙ（上記ではｙ∈Ｎ（ｘ）とされた。）との適合度を定量化し（上記では、ｐ（ｘ，ｙ）に基づいていた。）、該適合度が大きい場合には、画素値ｖ（ｙ）を利用した重み付き平均化処理において、当該画素ｙについて大きな寄与を認め、適合度が小さい場合には小さな寄与しか認めないようにすることで、当該画素ｘのノイズを有効に抑制する画像処理方法である、といえる。いわば、画素ｘと画素ｙとが「似たもの同士」のときには、該画素ｙを前記平均化処理に、より貢献させ、「似ていないもの同士」のときには、該画素ｙを殆ど又は全く無視する（重みをゼロ又はその近似値）、と言い換えてもよい。
【００５４】
このような処理を画像全体に施すことにより、画像のぼけ、つまり空間的分解能の低下を殆ど生じることなく、極めて高いノイズ抑制効果を発揮することができる。また、ノイズ抑制という用途に限定せず、例えばパターン認識の分野においても、重み関数、あるいはコヒーレントフィルタを好適な具体的形式にすることによって、優れた効果を発揮することができる。
【００５５】
ここで上記した「ダイナミックＣＴ」撮影とは、上記Ｘ線管１０１及びＸ線検出器１０２が被検体Ｐの同一部位を反復撮影（反復スキャン、連続回転型ＣＴ装置では、連続回転による反復撮影がしばしば行われる。）して次々に投影データを取得するとともに、該投影データに基づいて次々に再構成処理を行って時系列的な一連の画像を得る撮影方式のことをいう（この場合、画像表示部１０７における画像表示は、例えば図示しないカウンタ等によって、その画像の元となった投影データ収集に係るスキャン開始点又は終点から一定時間後に行われるように制御される。）。
【００５６】
したがって、このように取得・表示される画像は、映画等と同様に時系列的な複数枚の静止画像からなる、いわゆる動画像となる。なお、このような撮影方式は、典型的には、被検体Ｐに対し造影剤を注入し、その経時変化を観察・解析して、例えば血管における狭窄や閉塞等その他病変部の病態を分析するために用いられる。また、造影剤投与の前後２回だけに限り同一部位のＣＴ撮影を行う方式も、広義のダイナミックＣＴ撮影と考えることができる。
【００５７】
さて、従来においては、上記のような「ダイナミックＣＴ」撮影時、例えばＫ回の撮影を実施する間に被検体Ｐに何らかの変化（例えば、造影剤の濃度変化や呼吸動等が一般的に考えられる）があった場合、空間解像度を損なわず画像ノイズを抑制するためには、時間方向の平滑化を行うほかなかった。その結果、時間分解能が損なわれるという弊害は避け得なかった。
【００５８】
ところが、ダイナミックＣＴ撮影により取得される画像は、上述したように、動画像であって時間的変化を仔細に観察する目的で行うものであるから、その分解能が損なわれるというのは、本来、好ましい状況とは言えない。
【００５９】
コヒーレントフィルタを利用すれば、分解能を損ねず、Ｋ枚の静止画像のすべて（複数枚の画像）につきそのノイズを抑制することが可能な、次のようなダイナミック・コヒーレントフィルタ処理を実施することができる。
【００６０】
まず、上記のようにして得られた動画像たるＫ枚の静止画像につき定義される画素ｘについては、既に述べたように、画素値ｖ（ｘ）として、  
ｖ（ｘ）＝（ｖ１（ｘ），ｖ２（ｘ），…，ｖＫ（ｘ））  …  （１再掲）  
を構成することができる。ここで右辺各項における添え字１，２，…，Ｋは、Ｋ枚の各静止画像の通し番号である。
【００６１】
次に、この場合における重み関数ｗ１の具体的形式を、例えば次の（４）式により与える。
【数２】
2008100121.tif 000004
【００６２】
ただし、ｙ∈Ｎ（ｘ）であって、かつ、この集合Ｎ（ｘ）は、画素ｘにつき任意に設定してよい（＝どのような基準によって設定してもよい。）。しかし実際上は、画素ｘと該画素ｘから遠く離れた位置にある画素ｙとが仮説「ｖ（ｘ）＝ｖ（ｙ）。ただし、両画素のノイズに起因する差異を除く」を満たす可能性は一般に低いといえるから、集合Ｎ（ｘ）をｘに近接している画素の集合という基準で限定することは、演算速度向上等の実用的な意義がある。
【００６３】
したがってここでは、その一例として、集合Ｎ（ｘ）を、当該画素ｘを中心としたその周囲の矩形状エリアに含まれる画素の集合、とする。より具体的に、集合Ｎ（ｘ）としては、例えば、いま注目している静止画像一枚を構成する全画素が１２８×１２８画素であるような場合に、前記画素ｘを中心とした３×３画素分のエリアとしたり、また、５１２×５１２画素であるような場合に、当該画素ｘを中心とした１３×１３画素分のエリア等としてもよい。
【００６４】
また、上記（４）式におけるσｋは、ｋ枚目の静止画像の各画素が、そのどれにも共通な一定の程度で有するものと仮定して推定されたノイズの標準偏差であり、一方Ｃは、重みｗ１（ｐ（ｘ，ｙ））が、上記（４）式に代入された場合における作用の程度を決定調節可能なするパラメータである。
【００６５】
以下、これらσｋ及びＣについての説明を順に行う。  
まず、（４）式におけるσｋについて説明する（以下では、分散σｋ２として説明する。）。このσｋ２は、上述したように、ｋ枚目の静止画像上の各画素のスカラー値が有するノイズ成分の分散である。そしてまた、上記（４）式における分散σｋ２は、ｋ枚目の画像の各画素のスカラー値について一定値たる分散σｋ２を持つノイズを含んでいるものと仮定して推定したものである。一般に、このような仮定は、次に記すようなことを背景として、十分な正当性を持つ。
【００６６】
被検体Ｐの大きさ、Ｘ線管１０１及びＸ線検出器１０２、再構成部１０６等の構造が一定で、かつ、照射Ｘ線のエネルギを一定にした状態では、ＣＴ画像のノイズは、照射Ｘ線量、すなわちこれと比例関係にあるＸ線管１０１における管電流と照射時間との積（いわゆる管電流時間積（ｍＡ・ｓ））によって決定される。
【００６７】
一方、ＣＴ画像のノイズは加法的であり、概ねガウス分布に従うことも知られている。すなわち、ある画素ｘの画素値ｖ（ｘ）を構成する任意のスカラー値ｖｎ（ｘ）（ｎ＝１，２，…，Ｋ）について、その真値（ノイズの寄与分を除去した値）をｖｎ０（ｘ）とすると、これらの差の値ｖｎ（ｘ）−ｖｎ０（ｘ）は、概ね平均０、分散σｋ２のガウス分布に従う（なお、照射Ｘ線量ないし管電流時間積ｍ・Ａｓとノイズの分散σｋ２とは、概ね反比例関係にある。）。
【００６８】
また、この分散σｋ２は、画素ｘの位置そのもの（上で述べたように、例えば各座標値ｘ＝（ｘ，ｙ））にも依存するが、通常のＸ線ＣＴ装置１００においては、Ｘ線管１０１及びＸ線検出器１０２の間に、Ｘ線照射量を調節する物理的なＸ線フィルタ（例えば銅箔や金属塊等により構成された、いわゆる「ウェッジ」あるいは「Ｘ線フィルタ」と呼称されるもの）を備えているため、これを無視することができる。なぜならばウェッジは、被検体Ｐが水とほぼ同じ密度を持つ物質から構成されていることを利用して、どのＸ線検出器１０２においても同程度のＸ線量が検出されるよう、照射されるＸ線量の一部減弱する作用を有するものであり、従ってこのようなウェッジによれば、結果的に、ノイズの分散σｋ２を画素ｘの位置に殆ど依らない概ね一定値にする効果を生じるからである（ちなみに、このウェッジは、一般に、Ｘ線検出器１０２のダイナミックレンジを有効に利用することを本来の目的として設置されるものである。）。
【００６９】
以上のことから、ダイナミックＣＴ撮影により取得されたＫ枚の静止画像上においては、ｋ枚目の静止画像上におけるすべての画素について、分散σｋ２がほぼ一定であると推定することは妥当である。むろん、画素ごとに分散が異なる場合について本実施形態を拡張することも容易に推考できる。
【００７０】
さて次に、上記（２）式を具体的に演算するためには、その分散σｋ２として、どのような数値をあてるか、が問題となる。このようなことが問題となるのは、通常、ノイズの分布の形は想定できても（上記ではガウス分布）、分散σｋ２の具体値は不明であることが多いからである。
【００７１】
更に、一般的に、毎回の撮影毎に照射線量（Ｘ線管電流×照射時間（ｍＡｓ））を変更して撮影を行ってもよい。
【００７２】
さて、ｋ枚目の画像（ｋ＝１，２，…，Ｋ）に於いて各画素のスカラー値が持つノイズの分散をσｋ２とし、ｋ枚目の画像の撮影に用いた照射線量をＲｋとするとき、σｋ２はＲｋに比例する。従って少なくともひとつのｋ＝ｋ０についてσｋO２が指定できれば、他のｋに関しても、
【数３】
2008100121.tif 000005
【００７３】
によってσｋ２を正確に推定することができる。
【００７４】
本実施形態（このような事情が当てはまる）に於いては少なくともひとつのｋについて、以下のような方法でσｋ２の具体的数値の推定を行うことができる。
【００７５】
Ｋ回の撮影のうち、被検体Ｐに殆ど変化がなかったと仮定することのできるＮ回（１＜Ｎ≦Ｋ）の画像を用いて、実測により、分散σｋ２に対する期待値Ｅ[σｋ２]を求める方法が有効である。以下説明を簡単にするために、これらＮ枚の画像における照射線量は同じであり、従ってｋ＝１，２，…Ｎに関してσｋ２は一定（σ２と書く）と仮定する。これらＮ枚の画像における、ある画素ｘｆの画素値ｖ（ｘｆ）を構成する各スカラー値ｖ１（ｘｆ），ｖ２（ｘｆ），…，ｖＫ（ｘｆ）が含むノイズは、上述したように平均０、分散σ２のガウス分布に従うと予想されるから、これらの平均値を以下の（６）式、
【数４】
2008100121.tif 000006
【００７６】
を用いると、真の分散σ２に対する期待値Ｅ[σ２]を、
【数５】
2008100121.tif 000007
【００７７】
として求めることができる。そして、この分散の期待値Ｅ[σ２]は、上述した通り、Ｋ枚すべての静止画像上の全画素ｘにつき妥当するものと考えることができ、真の分散σ２の代用として用いるのに、一定程度以上確からしさが保証された値である。したがって、上記（４）式の実際の演算においては、このＥ[σ２]を（４）式のσ２に代入すればよい。
【００７８】
なお、このようなＥ[σ２]は、より具体的には、Ｋ枚の静止画像中、例えば１枚目と２枚目の静止画像に基づく実測値により求めてもよい（上記（６）及び（７）式で言えば、Ｎ＝２とすることに該当する。）。また、上記（６）及び（７）式の実際の演算に供される画素ｘｆについては、例えば、空気や骨が撮像されている部分を除いた適当な画素ｘｆのみを選定する（複数選定した場合は得られるＥ[σ２]すべての平均をとる）等といった工夫を施してもよい。さらに、その他一般的には、被検体Ｐの動きによる影響を抑える工夫等を施すと尚よい。
【００７９】
これらＮ枚の画像の撮影において照射線量が一定でない場合においても、σｋ２がＲｋに比例することを利用して正しくσｋ２を推定することは容易に推考できるであろう。
【００８０】
さて次に、上記（４）式におけるパラメータＣについての説明を行う。まず、（４）式においては、上記一般的形態で述べた危険率ｐ（ｘ，ｙ）の考え方が、以下のようにして含まれている。すなわち、（４）式の右辺分子における根号内の表式は、いわゆるχ二乗分布に従うとされる当該χ２値に一致するものであり、これを（２σ）２で除し、括弧の全体をｅの肩に置いた値は、危険率ｐ１（ｘ，ｙ）そのものである。つまり、
【数６】
2008100121.tif 000008
【００８１】
そして、上記（４）式は、この（８）式のように表されるｐ１（ｘ，ｙ）に関し、
【数７】
2008100121.tif 000009
【００８２】
としたものに他ならない。尚、Ａは定数でｐ１が（０〜１）の値になるように規格化されたものである。
【００８３】
結局、（４）式においては、上記したような一般的形態で述べた危険率ｐ（ｘ，ｙ）が陽には表示されてはいないが、重みｗ１（ｐ（ｘ，ｙ））の実態は、上述したように、まさしく危険率（＝ｐ１（ｘ，ｙ））の関数であると見ることができ（（９）式）、すなわち「適合度の関数」である（ただし、危険率と適合度とは、上述したように、一方が増えれば他方も増加する関係にある）。
【００８４】
そして、上記（９）式からわかるように、パラメータＣは、重みｗ１（ｐ（ｘ，ｙ））が、危険率ｐ１（ｘ，ｙ）にどの程度敏感に反応するかを決める効果がある。つまり、Ｃを大きくすると、ｐ１（ｘ，ｙ）がわずかに小さくなるだけで、ｗ１（ｐ（ｘ，ｙ））は０に近づく。また、Ｃを小さくするとそのような過敏な反応を抑制することができる。なお、Ｃとして、具体的には１乃至１０程度とすればよく、好適にはＣ＝３とするとよい。
【００８５】
本実施形態においては、両画素ｘ及びｙに関する類似判定、言い換えると、両画素ｘ及びｙに関する上述した帰無仮説Ｈの棄却の判定は、上述したことから明らかなように、上記危険率ｐ１（ｘ，ｙ）に基づいて、いわゆるχ二乗検定法（統計的検定法）によって決定されている。
【００８６】
また、上記（４）式の表式からわかるように、本発明においては、危険率ｐ（ｘ，ｙ）をｘ，ｙの組み合わせそれぞれについて計算した後、重みｗ（ｐ（ｘ，ｙ））を求めるといった手順を踏む必要は必ずしもなく、危険率ｐ（ｘ，ｙ）を具体的に求めずに、合成関数としての（ｗＯｐ）を、直接計算する構成としてもよい。
【００８７】
以上述べたように、分散σ２の推定をし（例えば、（７）式のＥ[σ２]）、かつ、パラメータＣを適当に決める（例えば、Ｃ＝３）ことにより、（４）式を用いて、ある画素ｘにつき定義される集合Ｎ（ｘ）（上述したように、例えば画素ｘを中心とした３×３画素分のエリア等）に含まれるすべての画素ｙについて、具体的な重みｗ１（ｐ（ｘ，ｙ））を求めることができる。後は、上記（２）式におけるｗ（ｐ（ｘ，ｙ））に代えて、このｗ１（ｐ（ｘ，ｙ））を用いることにより、コヒーレントフィルタの具体的な数値演算を実施することが可能となる。そしてその結果、時間分解能は勿論のこと、空間分解能をも損なわずに、ノイズを強く抑制した画素値ｖ´（ｘ）＝（ｖ´１（ｘ），ｖ´２（ｘ），…，ｖ´Ｋ（ｘ））（＝（３）式）、すなわちそのようなＫ枚の静止画像ないし動画像を、得ることができる。
【００８８】
このような画像処理を、概念的に把握しやすいよう図示したものが、図３である。すなわちまず、図３（ａ）においては、１，２，…，Ｋ枚ある静止画像において、ある画素ｘにつき、該画素ｘを中心とした３×３画素分の矩形状エリアＮ３×３（ｘ）が想定されている。この矩形状エリアＮ３×３（ｘ）の左角隅における画素を、ｙ１とすれば、この画素ｙ１は、図３に併せて示すように、画素値ｖ（ｙ１）を有している。
【００８９】
そして、この画素値ｖ（ｙ１）を構成するスカラー値ｖ１（ｙ１），ｖ２（ｙ１），…，ｖＫ（ｙ１）と画素値ｖ（ｘ）におけるスカラー値ｖ１（ｘ），ｖ２（ｘ），…，ｖＫ（ｘ）とのそれぞれにより、上記（４）式によって重みｗ１（ｐ（ｘ，ｙ１））が計算される（図３（ｂ））。また、矩形状エリアＮ３×３（ｘ）の残る画素ｙ２，…，ｙ８についても同様で、結局図３（ｂ）に示すように、ｗ１（ｐ（ｘ，ｙ１）），…，ｗ１（ｐ（ｘ，ｙ８））及び、ｗ１（ｐ（ｘ，ｘ））が得られる。（この場合、（８）式より危険率ｐ（ｘ，ｘ）は、“１”であり、したがって重みｗ１（ｐ（ｘ，ｘ））も、（９）式より“１”である（＝最大の重み付けがされている））。
【００９０】
次に、このようにして得られた、重みｗ１（ｐ（ｘ，ｙ１）），…，ｗ１（ｐ（ｘ，ｙ８）），ｗ１（ｐ（ｘ，ｘ））を、対応する画素の、ｋ枚目の画像におけるスカラー値ｖｋ（ｙ１），ｖｋ（ｙ２），…，ｖｋ（ｙ８），ｖｋ（ｘ）にそれぞれ乗算して総和を取り（上記（２）式における分子に該当する。）、これを矩形状エリアＮ３×３（ｘ）に関する重みｗ１の総和（同じく（２）式の分母に該当する。）により除せば、当該ｋ枚目の画像における画素ｘについての、ノイズが抑制されたスカラー値ｖ´ｋ（ｘ）を求めることができる（図３（ｃ））。また、ｋ＝１，２，…，Ｋのすべての画像につき、同じ重みｗ１（ｐ（ｘ，ｙ１）），…，ｗ１（ｐ（ｘ，ｙ８）），ｗ１（ｐ（ｘ，ｘ））を用いて、ノイズが抑制されたスカラー値ｖ´ｋ（ｘ）を求めることによって、画素ｘにおけるノイズが抑制された画素値ｖ´ｋ（ｘ）＝（ｖ´１（ｘ），ｖ´２（ｘ），…，ｖ´Ｋ（ｘ））が得られる。すべての画素ｘにつき、上記演算を繰り返せば、ノイズを抑制したＫ枚の画像が得られる。
【００９１】
このようにしてコヒーレントフィルタで算出された画素値ｖ´（ｘ）で構成される画像では、オリジナル画像で見られたランダムなノイズが、十分に抑制される。
【００９２】
なお、以上までに述べた各処理は、例えば図４に示すようなフローチャートに則ってこれを行えばよく、また、当該各処理に係る演算・画像表示等を実際のＸ線ＣＴ装置１００上で実現するためには、例えば、図２に示すように、分散値推定部１１１、重み演算部１１２及び画素値演算部１１３により構成される画像処理部１１０を設けて、これを実施すればよい。
【００９３】
このうち重み演算部１１２は、上述した手順通り、画素値ｖ（ｘ）及びｖ（ｙ）から直接重みｗ１（ｐ（ｘ，ｙ））を求める構成となっている。したがって当該演算部１１２は、危険率ｐ１（ｘ，ｙ）の値を具体的に求めることなく（すなわち、「危険率演算部（本発明にいう「適合度定量化部」を内蔵し）、重みを直接に求める装置である。なお、上記したような構成ではなく、具体的に危険率ｐ１（ｘ，ｙ）の値を求める「危険率演算部（適合度定量化部）」と、その出力に基づいて重みｗ１（ｐ（ｘ，ｙ））を求める「重み演算部」という、二段の手順を踏む構成としてもよい。いずれにせよ、重み演算部１１２は、分散値推定部１１１により推定された分散σ２と、ｖ（ｘ）及びｖ（ｙ）を用いて重みｗ１（ｐ（ｘ，ｙ））を算出する。
【００９４】
また、画素値演算部１１３は、画素値ｖ（ｘ）及びｖ（ｙ）、並びに重み演算部１１２により数値演算された重みｗ１（ｐ（ｘ，ｙ））を使って、画素値ｖ´（ｘ）を演算する。すなわち当該演算部１１３は、元となる画像のノイズを抑制する処理、すなわちコヒーレントフィルタの適用を実際に行う（以下、これを「コヒーレントフィルタをかける」と表現する。）。
【００９５】
上記のようなダイナミック・コヒーレントフィルタ処理においてＫ枚の静止画像から構成される動画像に、コヒーレントフィルタをかける場合には、上記画像処理部１１０における処理は、一旦すべての静止画像を再構成した後、これらを上記記憶装置１０Ｍに蓄え、後処理として後にこれらに対してコヒーレントフィルタをかけるようにしてもよいが、本実施形態はこのような形態に限定されるものではなく、上述した連続スキャン、連続投影データ収集、連続再構成及び連続表示という流れの中で、コヒーレントフィルタをかける処理をリアルタイムに実施する（以下、これを「リアルタイム・コヒーレントフィルタ処理」と呼ぶ。）のでもよい。
【００９６】
リアルタイム・コヒーレントフィルタ処理の好ましい実施形態においては、新しい画像が撮影され再構成されるたびに、以下のような処理を行う。最初に得られた画像（画像番号１）から最新の画像（画像番号Ｍ）までのうち、画像番号Ｍ，Ｍ−１，…，Ｍ−Ｋ＋１を持つＫ枚の静止画像上、共通の同一点（同一座標）の画素ｘの持つ画素値（スカラー値）を並べてＫ次元ベクトル値ｖ（ｘ）＝（ｖＭ（ｘ），ｖＭ−１（ｘ），…，ｖＭ−Ｋ＋１（ｘ））を構成する。こうして、上記の「ダイナミック・コヒーレントフィルタ処理」と全く同様にコヒーレントフィルタをかけることができる。ただし、画素値演算部１１３は実際には画素値ｖ´（ｘ）の全ての要素を計算するのではなく、最新の画像（画像番号Ｍ）に対応するスカラー値ｖＭ´（ｘ）だけを計算する。この結果、計算速度が向上するので、リアルタイムでノイズが抑制された最新の画像を表示できる。
【００９７】
この「リアルタイム・コヒーレントフィルタ処理」の別の好ましい実施例として、最初のＫ枚の画像が得られた時点で、上記と全く同様にコヒーレントフィルタをかけてｖ１´（ｘ），…，ｖＫ´（ｘ）を求めておき、以後は、Ｋ次元ベクトル値を画像番号Ｍ，Ｍ−１，…，Ｍ−Ｋ＋１を持つＫ枚の静止画像を用いてｖ（ｘ）＝（ｖＭ（ｘ），ｖＭ−１´（ｘ），…，ｖＭ−Ｋ＋１´（ｘ））によって構成し、これに対して上記のリアルタイム・コヒーレントフィルタ処理を適用するように構成してもよい。なお、これらのリアルタイム・コヒーレントフィルタ処理の際に画素値ベクトルｖ（ｘ）の次元Ｋを、マニュアル設定、あるいは自動設定によって、随時変更できるように構成しておくと便利である。
【００９８】
このようにコヒーレントフィルタにより、空間及び時間分解能を低下させることなく、ノイズだけを効果的に抑制したＣＴ画像を使ってＣＢＰスタディを実行し、脳組織（毛細血管）の血液の動態を定量的に解析し、そのインデックス（ＣＢＰ、ＣＢＶ、ＭＴＴ、Ｅｒｒ）を求めることにより、その精度及び信頼性の向上が期待できる。
【００９９】
以下に、ＣＢＰスタディについて説明する。  
（ＣＢＰスタディ）  
（原理）  
ＣＢＰスタディでは、脳組織内の”毛細血管を通過する血流”の動態を定量的に表すＣＢＰ、ＣＢＶ、ＭＴＴ、Ｅｒｒのインデックスを求め、またこれらインデックスのマップを出力する。
【０１００】
ＣＢＰ：脳組織の毛細血管内の単位体積及び単位時間あたりの血流量［ml/１００ml/miｎ］  
ＣＢＶ：脳組織内の単位体積あたりの血液量［ml/１００ml］  
ＭＴＴ：毛細血管の血液平均通過時間［秒］  
Ｅｒｒ：解析モデルからの実測値のずれ残差の指標。この指標の多少により、脳動脈の支配組織と非支配組織との識別等の解析が可能となる。
【０１０１】
ＣＢＰスタディでは、トレーサーとして脳血管透過性を持たない造影剤、たとえばヨード造影剤が使用される。ヨード造影剤は例えばインジェクターにより肘静脈から注入される。インジェクターにより静注されたヨード造影剤は、心臓、肺を経由して、脳動脈から流れ込む。そして、脳動脈から、脳組織内の毛細血管を経て、脳静脈へと流れ出ていく。このとき、脳血管透過性を持たない造影剤、たとえばヨード造影剤は正常な脳組織内の毛細血管では造影剤は血管外へ漏れ出ることなく通過する。
【０１０２】
造影剤の通過の様子をダイナミックＣＴで連続的に撮影して、その連続画像から、脳動脈上の画素の時間濃度曲線Ｃａ（ｔ）、毛細血管を含む脳組織上の画素の時間濃度曲線Ｃｉ（ｔ）、脳静脈上の画素の時間濃度曲線Ｃｓｓｓ（ｔ）をそれぞれ測定する。
【０１０３】
ＣＢＰスタディでは、造影剤の血中濃度について脳組織に近い脳血管の血中濃度の時間曲線Ｃａ（ｔ）と、毛細血管の血中濃度の時間曲線Ｃｉ（ｔ）との間で成り立つ理想的な関係を解析モデルとしており、つまり脳組織に入る直前の血管から造影剤を注入した場合、毛細血管を含む脳組織単位体積（１画素）内の時間濃度曲線は立ち上がりが垂直で、若干の勾配を持って立ち下がる形になる。これは、矩形関数で近似することができる（ｂｏｘ−ＭＴＦ法：ｂｏｘ−Ｍｏｄｕｌａｔｉｏｎ  Ｔｒａｎｓｆｅｒ  Ｆｕｎｃｔｉｏｎ  ｍｅｔｈｏｄ）。
【０１０４】
脳動脈血中時間濃度曲線Ｃａ（ｔ）を入力関数、脳組織の時間濃度曲線Ｃｉ（ｔ）を出力関数として、毛細血管を通過する過程の特徴を、矩形関数で表される伝達関数として求めることができる。
【０１０５】
（具体的な手順）
図５、図６には、本実施形態によるＣＢＰスタディの典型的手順を示している。  
まず、肘静脈等の血管にボーラスインジェクション（造影剤を一気に投与する）を行い、その直後あるいは直前からダイナミックＣＴ（同じ箇所を反復して撮影する）を行う。最も典型的な手技として、肘静脈へボーラスインジェクションを行った場合、概ね２０〜４０間に渡って、例えば０．５〜２秒間隔で撮影を繰り返す。ダイナミックＣＴで得たＮ枚のＣＴ画像のうちのｊ枚目の各ピクセル（ｘ、ｙ）のＣＴ値をｖ（ｘ、ｙ、ｊ）とする。これはこの画素（ｘ、ｙ）における時間濃度曲線（滑らかな曲線である）ｆ（ｔ、ｘ、ｙ）をサンプリングしたものに他ならない。
【０１０６】
まず、前処理として、ステップＳ１で、ＣＴ画像各々から、明らかに脳組織以外の組織であることが判別される画素を、解析対象から除外する。すなわち、脳組織のＣＴ値として考えられる範囲（例えばＣＴ値１０〜６０ＨＵ）に入らない値を示す画素は、空気や骨、脂肪などに対応する画素であり、脳血流の定量とは関係ないのでこれらは無視して良い。この解析範囲は、デフォルトとして、１０〜６０ＨＵに設定されるが、入力部１０９を介して任意に設定可能である。
【０１０７】
また、前処理として、ステップＳ２で、造影効果の初期化が行われる。各画素に於ける造影効果（ＣＴ値の上昇）を得るためには、各画素（ｘ，ｙ）について、その画素に対応する組織に造影剤が到達する以前の画像（一般に複数枚得られる）を、通し番号１，２，…Ｋで表すと、その時間的平均値は、
【数８】
2008100121.tif 000010
【０１０８】
を求め、この値をｂ（ｘ，ｙ）とする。そして、ｊ＝Ｋ＋１、Ｋ＋２、…，Ｎの各画像の画素値ｖ（ｘ，ｙ，ｊ）について、  
ｑ（ｘ，ｙ，ｊ）＝ｖ（ｘ，ｙ，ｊ）−ｂ（ｘ，ｙ）  
ｊ＜Ｋについて  
ｑ（ｘ，ｙ，ｊ）＝０  
とすればよい。処理を簡単にするためには、どの画素に関しても同じＫを採用しても良い。こうして得られたｑ（ｘ，ｙ，ｊ）は、造影効果の（滑らかな）時間濃度曲線ｑ（ｔ，ｘ，ｙ）をｔ＝ｔ１，ｔ２，…ｔＮにおいてサンプリングしたものに他ならないと考えることができる。このｑ（ｔ，ｘ，ｙ）を用いて脳血流の定量解析を行う。
【０１０９】
定量解析にあたってはまず、右脳エリアと左脳エリアをＣＴ画像上で分離する必要がある。上述したようにＣＢＰスタディでは、毛細血管の血流動態の様子を、脳動脈の時間濃度曲線Ｃａ（ｔ）に対する脳組織の時間濃度曲線Ｃｉ（ｔ）の伝達関数ＭＴＦとして求めるものであり、従って、解析対象の脳組織が、Ｃａ（ｔ）の脳動脈の支配下にがあることが前提とされる。少なくとも左脳と右脳とでそれぞれ別々の脳動脈の時間濃度曲線Ｃａ（ｔ）を使って個別に解析する、つまり左脳の脳動脈の時間濃度曲線Ｃａ（ｔ）は同じ左脳の脳組織の解析にだけ使用し、動揺に、右脳の脳動脈の時間濃度曲線Ｃａ（ｔ）は同じ右脳の脳組織の解析にだけ使用することが必要とされる。
【０１１０】
脳を左脳エリアと右脳エリアとを分割するために、図７に示すように、ＣＴ画像上に分割線が画面上に図形として重ねて表示される（Ｓ３）。分割線が最初は画像中央に表示されるように構成しても良い。操作者は、画像を参照して、分割線を移動し、また分割線を構成する複数の構成点を移動して任意に屈曲させることにより、左右エリアを分割する。
【０１１１】
次に、ＣＴ画像上で脳動脈上に脳動脈ＲＯＩを設定するのであるが、この設定の精度を向上し且つ容易にするためにＲＯＩ設定支援部１２１により支援マップが作成され、ＣＴ画像とは別に、またはＣＴ画像に重ねて表示される（Ｓ４）。支援マップとしては、例えば、ＡＴ（アピアランスタイム）マップ、ＰＴ（ピークタイム）マップ、ＴＴ（トランジットタイム）マップがあげられる。各画素について、図８に示すように、造影前の任意の時刻（例えばデータ収集開始時刻）Ｔ０から造影剤濃度がピークｐｅａｋの数パーセント（例えば１パーセント）に達する時刻までの時間ＡＴ、時刻Ｔ０から造影剤濃度がピークに達した時刻までの時間（ピークタイム）ＰＴ、又は造影剤の移動時間を例えば半値幅で表すＴＴが計算され、マップとして生成され表示される。デフォルトでは、これらＡＴマップ、ＰＴマップ、ＴＴマップの全種類が生成され表示されるようになっているが、任意の１種類、又は任意の２種類を操作者が選択することが可能である。
【０１１２】
脳動脈のこれら数値は他と比較して高いある程度固有値で現れる傾向にあるので、その値を中心としたウインドウ内が抽出して表示されるよう設定されたカラーマップ（ルックアップテーブル）を通してカラー表示させることで、脳動脈の場所を容易に識別し、脳動脈ＲＯＩを正確に設定することが可能となる（Ｓ５）。典型的には、左右脳エリアそれぞれに、脳動脈ＲＯＩは、前大脳動脈（ＡＣＡ）、中大脳動脈（ＭＣＡ）、後大脳動脈（ＰＣＡ）の３箇所ずつ設定される。
【０１１３】
なお、マルチスライス、撮影の場合、例えば隣接する４枚のスライスを解析対象とする場合、図１２に示すように、その各スライスで個々に脳動脈ＲＯＩを設定することは作業負担が大きいばかりで、解析を行う上では必要のない作業である。従って、ある任意の１スライスで設定した脳動脈ＲＯＩを他のスライスにも共用する。または、後述するコヒーレントレグレッション法を使って全スライスで共通に用いることのできる脳動脈の時間濃度曲線Ｃａ（ｔ）を作成するようにしてもよい。
【０１１４】
次に、設定された脳動脈ＲＯＩ各々に関して時間濃度曲線Ｃａ（ｔ）が、ダイナミックＣＴによる連続画像データから時間濃度曲線作成部１２２により作成される（Ｓ６）。
【０１１５】
ここで、画素に比べて、脳動脈は非常に細いものが多く、しかも一般にＣＴの撮影スライスに対して直交していないため、画像上のどの画素も正確に動脈血のＣＴ値を表しておらず、１画素が脳動脈と他の組織との混在により構成され、そのパーシャルボリューム効果のためにそれよりも低い造影効果しか示さないことがほとんどである。また、動脈をパーシャルボリュームとして含むこれらの画素において、画像ノイズが大きい。特に脳梗塞を生じている部位の動脈等においては、造影効果が比較的小さいために、ノイズの影響は甚大である。画像ノイズに関しては上述したコヒーレントフィルタにより抑制されているものの、パーシャルボリューム効果の影響は依然として残存している。
【０１１６】
この問題は、脳動脈の時間濃度曲線は、単一のスライス画像において計測するのではなく、その動脈を含む立体内の画素を用いて後述するコヒーレントレグレッション法を適用することで抑圧することが可能である。従って、上述したコヒーレントフィルタ法に代えて、この段階でコヒーレントレグレッション法を適用するようにしてもよい。
【０１１７】
また、この方式によれば、動脈ごとにそれに対応する唯一の脳動脈の時間濃度曲線が得られ、従って、撮影範囲内にある全スライス中の任意の部位の解析に利用することができ、これによって、特定の動脈について、その脳動脈の時間濃度曲線が最も明瞭に得られるスライスを選んで、その脳動脈の時間濃度曲線を全スライスに適用することができ、脳動脈の時間濃度曲線の数を減らすことができる。
【０１１８】
（コヒーレントレグレッション法）  
上記時間濃度曲線の作成では、パーシャルボリューム効果、ランダムノイズの影響を除去することが重要である。まず、「時間濃度曲線」とは、上記ダイナミックＣＴ画像中の特定の部位における像の濃度値の経時的変化を表す曲線である。ことに、上記医用画像診断装置においては、人体組織等における血流動態や代謝機能等の詳細を調べる事を目的として、人体の特定組織内の造影剤濃度等の経時的変化を時間濃度曲線として計測することが行われている。また、天体観測等においては、特定の天体の光度変化等を解析する目的で、時間濃度曲線が用いられる。より形式的に明示すると、すなわち、時間濃度曲線とは、時刻ｔｋにおけるある部位の濃度値をｄｋとするとき、対の列｛＜ｔｋ，ｄｋ＞（ｋ＝１，２，・・・，Ｋ）｝として表現される。また、時間濃度曲線の多くの用途においては、必ずしもｄｋの絶対的な値が必要なのではなく、むしろ最初の画像１を基準とする増分（ｄｋ−ｄ１）だけが得られれば十分である。さらにそのような用途のうちの多くでは、単に（ｄｋ−ｄ１）に比例するデータ Ａ（ｄｋ−ｄ１）（ここにＡは未知の比例係数）だけが得られれば十分である。この場合には、従って、対の列｛＜ｔｋ，Ａ（ｄｋ−ｄ１）＞（ｋ＝１，２，・・・，Ｋ）｝が、求める時間濃度曲線である。
【０１１９】
このような時間濃度曲線を求めるためには、原理的には、上記ダイナミックＣＴ画像を構成する各画像ｋ（ｋ＝１，２，・・・，Ｋ）における、該時間濃度曲線を測定しようとする部位に含まれる画素ｘのスカラー値ｖｋ（ｘ）を用いて、対の列｛＜ｔｋ，ｖｋ（ｘ）＞｝あるいは、｛＜ｔｋ，Ａ（ｖｋ（ｘ）−ｖ１（ｘ））＞｝を構成すればよい。
【０１２０】
しかし、実用においては、上記医用画像診断装置等によって撮影されたダイナミックＣＴ画像にランダムなノイズが含まれているために、本来測定しようとする時間濃度曲線を正確に求められないという問題がある。
【０１２１】
さらに、実用においては、これらのダイナミックＣＴ画像においては、いわゆる「パーシャルボリューム効果」が生じる。パーシャルボリューム効果とは、すなわち、被検体内の微小な物体の像は、画像上では少数個の画素によって表現されるが、これら少数個の画素には、被検体内の隣接する物体の像も影響を与えるため、これら少数個の画素の画素値は（本来計測しようとする濃度値の変動に比例するものの）比較的小さな変動しか示さない、という現象である。言い換えれば、これら少数個の画素の画素値は僅かな信号しか含まない。従って、パーシャルボリューム効果が生じている場合には、どの画素ｘを取っても対の列｛＜ｔｋ，ｖｋ（ｘ）＞（ｋ＝１，２，・・・，Ｋ）｝は非常に信号レベルが低く、本来計測しようとしているのではない組織における濃度値の変化の影響を受け、さらにランダムなノイズが存在するために、本来測定しようとする時間濃度曲線｛＜ｔｋ，ｄｋ＞｝を正確に求められないという問題がある。
【０１２２】
そこで、従来は、ランダムなノイズを抑制するために、時間的又は空間的平滑化を用いていたが、時間平均を行うと時間分解能が損なわれてしまい、また、空間平均を行うと、本来の該測定しようとする部位以外の部位の濃度の経時変化が計測値に混入するという問題点があった。このような問題点を解決し、より正確な時間濃度曲線を得るために、コヒーレントフィルタを採用する。
【０１２３】
まず、本実施形態のコヒーレントフィルタにおいて用いるべき、帰無仮説について説明する。計測しようとする部位における真の時間濃度曲線を｛＜ｔｋ，ｄｋ＞（ｋ＝１，２，・・・，Ｋ）｝であると仮定するとき、その一次変換である｛＜ｔｋ，Ａ（ｄｋ−ｄ１）＞（ｋ＝１，２，・・・，Ｋ）｝（だたしＡは未知の係数）を計測することを目的とする場合において、計測しようとする部位に概ね相当する画素の集合Ｒを設定する。この集合Ｒの要素である任意の画素ｘ∈Ｒについて、条件Ｑ：「もし、この画素ｘが上記真の時間濃度曲線を良く反映し、しかも他の部位の経時的濃度変化の影響をほとんど受けていない」のであれば、（ベクトル値としての）画素値ｖ（ｘ）＝（ｖ１（ｘ），ｖ２（ｘ），．．．，ｖＫ（ｘ））について、パーシャルボリューム効果およびランダムノイズの影響を考慮することによって、  
ｖｋ（ｘ）＝ｐ（ｘ） ｄｋ＋ｑ（ｘ）＋γｋ（ｘ）          …（１１）  
（ｋ＝１，２，・・・，Ｋ）  
が成り立つと仮定することができる。ここに、ｐ（ｘ）およびｑ（ｘ）は、画素ｘごとに異なるが画像番号ｋ（すなわち撮影時刻ｔｋ）によっては変化しない未知の係数であり、パーシャルボリューム効果をモデル化したものである。またγｋ（ｘ）はランダムなノイズをモデル化したものであって、画素ｘごとに、しかも画像番号ｋごとに値が異なるが、その期待値は０であり、またその統計分布は画素ｘにも画像番号ｋにも依存しない。
【０１２４】
以上の仮定によれば、該集合Ｒの要素である任意の２個の画素ｘ，ｙに関して、もし「画素ｘ，ｙが共に（上記の）条件Ｑを満たす。」という命題が成り立つのであれば、次式の関係が成り立つことが証明できる。
【０１２５】
ｖｋ（ｘ）＝ａ１ｖｋ（ｙ）＋ａ２＋ξｋ    （ｋ＝１，２，・・・，Ｋ）  
…  （１２）
ここに、ａ１およびａ２は、画素の組ｘ，ｙごとに異なるが画像番号ｋ（すなわち撮影時刻ｔｋ）によっては変化しない未知の係数である。またξｋはランダムなノイズであって、画素の組ｘ，ｙごとに、しかも画像番号ｋごとに値が異なるが、その期待値は０である。
【０１２６】
（１２）式は以下のようにして導かれる。すなわち、ｘにｙを代入して得られる式  
ｖｋ（ｙ）＝ ｐ（ｙ） ｄｋ＋ｑ（ｙ）＋γｋ（ｙ）…  （１３）
を変形すると、
【数９】
2008100121.tif 000011
【０１２７】
とおくことによって、（１２）式が導かれる。ここで、（１６）式のａ１とａ２は、パーシャルボリューム効果を現すパラメータであり、また（１６）式のξｋは、ランダムなノイズを表す。
【０１２８】
以上から、「画素ｘ，ｙが共に条件Ｑを満たす。」という命題は、帰無仮説Ｈ０´「ｖｋ（ｘ）＝ａ１ｖｋ（ｙ）＋ａ２＋ξｋ    （ｋ＝１，…，Ｋ）である。」と等価であることが示された。
【０１２９】
次に帰無仮説Ｈ０´「ｖｋ（ｘ）＝ａ１ｖｋ（ｙ）＋ａ２＋ξｋ    （ｋ＝１，…，Ｋ）である。」を、実質的に等価であり、かつ実際に検定できる形式の命題に変換する方法について述べる。この帰無仮説を改めて数学的に厳密な表現で述べると、帰無仮説Ｈ０´「ある定数ａ１およびａ２が存在して、ξｋ＝ｖｋ（ｘ）−ａ１ｖｋ（ｙ）−ａ２（ｋ＝１，…，Ｋ）は平均０、分散（σ２ｈ（ａ１））の正規分布に従う。」となる。ここに係数ｈ（ａ１）は、  
ｈ（ａ１）＝１＋ａ１２                            …  （１７）
である。（（１７）式はａ１とξｋの定義である（１６）式、および、ランダム変数に関する分散の持つ一般的な性質から直ちに導かれる。）また、上記の分散σ２の値は、簡便かつ実用上十分正確に推定できる。
【０１３０】
以上から、もし、上記の定数ａ１およびａ２を決定することができれば、上記の帰無仮説Ｈ０´を検定することが可能である。そして実際上は、これらの定数の最適な推定値ａ１〜およびａ２〜が得られれば十分である。
【０１３１】
このような、定数ａ１およびａ２の最適な推定値の算出には、公知の当てはめ法（fitting）がそのまま利用できる。そこで、以下では、そのような当てはめ法の典型的な具体例として、線形最小二乗法を用いる場合における概要を説明する。線形最小二乗法を本実施例に適用するには、単に、上記の帰無仮説のξｋの二乗和をＳ（ａ）として、すなわち
【数１０】
2008100121.tif 000012
【０１３２】
を定義する。Ｓ（ａ）の値は定数ベクトルａ＝（ａ１，ａ２）、すなわち上記の定数ａ１およびａ２、の値に依存する。このＳ（ａ）が最小の値を取るような定数ベクトルａを算出すれば、定数ａ１およびａ２に関する、不偏推定の意味での最適な推定値ａ１〜およびａ２〜が得られる。なお、線形最小二乗法の具体的な計算方法としては、様々な公知の方法を利用することができ、しかも、これら公知の計算方法はいずれも非常に簡単であり、必要な計算時間はごく僅かである。
【０１３３】
このようにして、上記の定数ａ１，ａ２の最適な推定値ａ１〜，ａ２〜を算出した結果、次式で定義される残差  
ｒｋ〜（ｘ，ｙ）＝ｖｋ（ｘ）−ａ１〜ｖｋ（ｙ）−ａ２〜        …（１９）
を具体的に計算することができる。従って、この残差ｒｋ〜を用いて、上記の帰無仮説Ｈ０´を、実質的に等価な帰無仮説Ｈ０”「ｒｋ〜（ｘ，ｙ）（ｋ＝１，…，Ｋ）は平均０、分散（１＋（ａ１〜）２）σ２の正規分布に従う。」と言い換えることができる。これは、実際に検定の計算を実行可能な具体的命題である。
【０１３４】
なお、さらに、ベクトルによる表現
【数１１】
2008100121.tif 000013
【０１３５】
（ただし、ベクトルａ及びξは画素の組ｘ，ｙに依存する。）を導入し、また、次式
ｆ（ａ~，ｖ（ｙ））＝ａ１〜ｖ（ｙ）＋ａ２〜   …  （２１）
で定義されるベクトル値関数ｆを用いて帰無仮説Ｈ０´を言い換えると、帰無仮説Ｈ０”は「ｖ（ｘ）＝ｆ（ａ~，ｖ（ｙ））＋ξ、（ただし、ξは平均０、分散（１＋（ａ１〜）２）σ２の正規分布に従う。）」となり、これは上述した帰無仮説Ｈ０と全く同じ形式である。すなわち、本実施形態は上述したコヒーレントフィルタの一変形例であることは明らかである。なお、ここで、上記ｆ（ａ~，ｖ（ｙ））とは、すなわち、画素ｙの画素値ｖ（ｙ）に対して、パーシャルボリューム効果を現すパラメータａを最適に調節して、画素ｘの画素値ｖ（ｘ）と最も高い適合度を持つように変換したものを意味する。
【０１３６】
次に、本実施形態において、上記の帰無仮説Ｈ０”を用いて、コヒーレントフィルタによって時間濃度曲線を求める方法について説明する。計測しようとする部位に概ね相当する画素の集合Ｒについて、この集合Ｒに含まれるあるひとつの画素ｘ∈Ｒについて、集合Ｒの要素である全ての画素y∈Ｒに対して、以下の計算を行う。すなわち、上記の方法を用いて実際に残差ｒｋ〜（ｘ，ｙ）（ｋ＝１，…，Ｋ）を算出し、次に、上記の帰無仮説Ｈ０”「ｒｋ〜（ｘ，ｙ）（ｋ＝１，…，Ｋ）は平均０、分散（１＋（ａ１〜）２）σ２の正規分布に従う。」を棄却する場合の危険率ｐ（ｘ，ｙ）ないし重みｗ（ｐ（ｘ，ｙ））を具体的に計算する。そして、重み付き平均ｖ´ｋ（ｘ）を下式（２２）によって計算し、画素ｘにおける時間濃度曲線｛＜ｔｋ，ｖ´ｋ（ｘ）−ｖ´１（ｘ）＞（ｋ＝１，２，・・・，Ｋ）｝を構成する。
【数１２】
2008100121.tif 000014
【０１３７】
こうして得られた時間濃度曲線は、画素ｘにおける真の時間濃度曲線｛＜ｔｋ，ｄｋ＞｝の一次変換である｛＜ｔｋ，Ａ（ｄｋ−ｄ１）＞｝（だたしＡは未知の係数）を近似している計測値であり、しかも、重み付き平均の効果によって、ランダムなノイズが抑制されている。また、他の画素ｙの画素値ベクトルに対しては、式から明らかなように、パーシャルボリューム効果の影響を補正したものが用いられている。さらに、本実施形態はコヒーレントフィルタの共通の特徴である、「時間平均を全く使用せず、また空間平均を画素ｘとの適合度に基づく重みを使って計算する」という性質を有する。従って、本実施形態によって、時間分解能を損なわず、パーシャルボリューム効果の影響を抑制し、しかもランダムなノイズが抑制された時間濃度曲線を得ることができる。なお、このようにして時間濃度曲線を求める方式を、特に「コヒーレントレグレッション法」と称す。
【０１３８】
次に、具体的に、医療用のＸ線ＣＴにおけるダイナミックＣＴ撮影等で得られたダイナミックＣＴ画像における、時間濃度曲線の臨床的利用の一例を説明する。この応用例では、造影剤を血管に急速に注入しながら、ダイナミックＣＴ等の撮影を行い、人体組織中に存在する動脈の像の濃度変化を時間濃度曲線として計測することによって、当該組織における血流動態を診断しようとするものである。
【０１３９】
この応用例において、多くの場合、人体組織中の動脈は一般に非常に細いために、ＣＴによる断層画像上に現れる動脈の像は、パーシャルボリューム効果を生じる。さらに、像にはランダムなノイズが含まれていることは言うまでもない。このため、従来の方法では、動脈に関する十分に正確な時間濃度曲線を得ることは困難であり、強いて計測を行えば、動脈に関する真の時間濃度曲線＜ｔｋ，Ｄｋ＞の一次変換である＜ｔｋ，Ａ（Ｄｋ−Ｄ１）＞（ここにＤｋは動脈の像に相当する一群の画素の、時刻ｔｋにおける（スカラー値である）画素値を表す。また、ｋ＝１，２，・・・，Ｋ）をある程度近似する測定値＜ｔｋ，（ｖｋ（ｘ）−ｖ１（ｘ））＞しか得られなかった。この測定値はランダムなノイズを含む。また、パーシャルボリューム効果の影響のために、係数Ａは未知のままである。
【０１４０】
そこで、＜ｔｋ，Ａ（Dｋ−D１）＞を十分に近似する測定値＜ｔｋ，（ｖ´ｋ（ｘ）−ｖ´１（ｘ））＞（ｋ＝１，２，・・・，Ｋ）を得ることができる。一方、同じ断層画像上で観察できる静脈の中には、相当に太いものが存在し、従ってそれらの静脈に関しては、従来の方法で、時間濃度曲線の十分に良い近似値＜ｔｋ，（Ｊｋ−Ｊ１）＞（ｋ＝１，２，・・・，Ｋ）を得ることができる。ここにＪｋは静脈の像に相当する一群の画素の、時刻ｔｋにおける画素値を表す。
【０１４１】
ところで、血液循環に関する時間濃度曲線においては、命題Ｓ:「もし、時刻ｔ１における血中の造影剤濃度が０であるならば、どの血管ｄに関する時間濃度曲線＜ｔｋ，（ｄｋ-ｄ１）＞も、その曲線下面積（ＡUC：Ａrea Under CurＶe）が一致する」という性質が成り立つことが知られている。ここで言う曲線下面積とは、時間濃度曲線＜ｔｋ，（ｄｋ-ｄ１）＞の時間tに関する積分を意味する。
【０１４２】
従って、ある血管ｄに関する時間濃度曲線＜ｔｋ，（ｄｋ-ｄ１）＞の曲線下面積ＡUC（ｄ）は、例えば次式によって近似的に計算することができる。
【数１３】
2008100121.tif 000015
【０１４３】
従って、静脈に関して従来の方法で得られた時間濃度曲線｛＜ｔｋ，（Ｊｋ-Ｊ１）＞｝に関する曲線下面積ＡUC（Ｊ）を（２２）式を用いて計算することができる。（ｄにＪを代入すればよい。）また、動脈に関して、仮に、時間濃度曲線｛＜ｔｋ，（Dｋ−D１）＞｝が知られていれば、曲線下面積ＡUC（D）を（１８）式を用いて同様に計算することができ、しかも上記命題Sに従って
ＡUC（D）≒ＡUC（Ｊ）                                       …  （２４）
が成り立つ筈である。しかし、実際には、時間濃度曲線＜ｔｋ，（Dｋ−D１）＞は未知であるため、ＡUC（D）は計算できない。
【０１４４】
一方、本発明に係る方式で得られた時間濃度曲線＜ｔｋ，（ｖ´ｋ（ｘ）−ｖ´１（ｘ））＞は、＜ｔｋ，Ａ（Dｋ−D１）＞を近似するものであり、後者は未知の係数Ａを含んでいる。このため、｛＜ｔｋ，（ｖ´ｋ（ｘ）−ｖ´１（ｘ））＞｝から（２３）式を用いて具体的に計算できる曲線下面積ＡUC（ｖ´）は、ＡUC（D）のちょうどＡ倍でなくてはならない。すなわち、  
ＡUC（ｖ´）≒Ａ ＡUC（D）                               …  （２５）  
である。すなわち、（２４）式と（２５）式から、  
Ａ≒ＡUC（ｖ´）／ＡUC（Ｊ）                            …  （２６）  
という関係が成り立つ。（２６）式の右辺は（２３）式を用いて具体的に計算できるため、未知であった係数Ａの値が具体的に決定できる。そこで、この係数Ａの値を用いて時間濃度曲線＜ｔｋ，（ｖ´ｋ（ｘ）−ｖ´１（ｘ））／Ａ＞を構成すれば、これは、動脈の時間濃度曲線＜ｔｋ，（Dｋ−D１）＞を近似するものに他ならない。このように、曲線下面積を用いて、未知であった比例係数Ａの値を決定した時間濃度曲線を構成する方法を「ＡUC法」と呼ぶ。
【０１４５】
以上から、ダイナミックＣＴ撮影等で得られたダイナミックＣＴ画像における、時間濃度曲線の臨床的利用において、上記コヒーレントレグレッション法に、さらに上記ＡUC法を組み合わせることによって、従来の方法では計測が困難あるいは不可能であった、細い動脈の時間濃度曲線に関しても、パーシャルボリューム効果およびランダムなノイズの影響を排除し、しかも、未知の比例係数Ａを含まない測定値が得られる。
【０１４６】
なお、もちろんＡUC法は、単独で従来の方法で計測された動脈に関する時間濃度曲線＜ｔｋ，（ｖ´ｋ（ｘ）−ｖ´１（ｘ））＞に対しても適用でき、（ランダムなノイズやパーシャルボリューム効果の影響は排除できないものの、）未知であった比例係数Ａの値を決定した時間濃度曲線を構成できる。
【０１４７】
（上矢状静脈洞の時間濃度曲線を使った脳動脈の時間濃度曲線の補正（パーシャルボリューム効果の影響を抑圧））  
パーシャルボリュームの影響を抑圧するために、このコヒーレントレグレッションに代えて又は併用して、上矢状静脈洞の時間濃度曲線Ｃｓｓｓ（ｔ）を使って脳動脈の時間濃度曲線Ｃａ（ｔ）を補正するようにしてもよい。
【０１４８】
まず、図５のステップＳ７において、図１１に示すように、ＣＴ画像上で上矢状静脈洞を取り囲むように大き目に上矢状静脈洞ＲＯＩを設定する。上矢状静脈洞は動脈に比べて大きく、また位置も固定的であるので、上矢状静脈洞ＲＯＩ設定は容易である。この大き目の上矢状静脈洞ＲＯＩには、複数の画素が含まれている。
【０１４９】
次に、上矢状静脈洞ＲＯＩの中の全ての画素がその全域に渡って上矢状静脈洞に含まれるように、上矢状静脈洞ＲＯＩを縮小処理する（Ｓ８）。縮小処理としては、例えば、まず、上矢状静脈洞ＲＯＩ内の画素各々について、しきい値処理（二値化）を実行し、ＲＯＩ内の二値マップ（“０”，“１”）を作成する。当該しきい値は、上矢状静脈洞の像を、その周辺の組織や骨の像と分離する値に設定される。“１”は上矢状静脈洞の像上の画素であることを表し、“０”は周辺の組織や骨の像上の画素であることを表している。この二値マップの各画素（中心画素）を、その近傍４または８画素の値に従って、置換する。中心画素が“１”であって、且つ、近傍４または８画素全てが“１”である場合のみ、中心画素の値を“１”に維持する。つまり、中心画素が“０”である場合はもちろん、たとえ“１”であったとしても、近傍４または８画素のなかの１つでも“０”を示しているときには、当該中心画素の値を“０”に置換する。従って、上矢状静脈洞ＲＯＩは、上矢状静脈洞の像の外形よりも、少なくとも１画素分縮小される。それにより縮小処理を受けた上矢状静脈洞ＲＯＩの中の全ての画素は、上矢状静脈洞像上の画素であるという条件を高い確度で実現され得る。
【０１５０】
また、この手法に代えて、時間濃度曲線の曲線下面積ＡＵＣを使って上矢状静脈洞ＲＯＩを修正するようにしてもよい。この場合、大き目のＲＯＩを探索範囲として、その中の画素各々について時間濃度曲線の曲線下面積ＡＵＣを計算する。造影効果により上矢状静脈洞像上の画素の曲線下面積ＡＵＣは、周辺画素のそれに比べて明らかに高位を示す。従って、この曲線下面積ＡＵＣに対してしきい値処理を実行することにより、ＲＯＩの中から上矢状静脈洞像上の画素だけを選別することができる。
【０１５１】
このようにしていずれかの手法又は両手法を併用してアンド条件によりピックアップされた上矢状静脈洞像上であるという確度の高い複数の画素に対して、各画素の時間濃度曲線が平均化され、上矢状静脈洞の時間濃度曲線Ｃｓｓｓ（ｔ）が作成される（Ｓ９）。
【０１５２】
ここで、ヨード造影剤は関門非通過であるので、原理的に、ヨード濃度は脳動脈と脳静脈とで変化しない、つまり、上矢状静脈洞の時間濃度曲線Ｃｓｓｓ（ｔ）の曲線下面積ＡＵＣは、Ｓ６で作成した脳動脈の時間濃度曲線Ｃａ（ｔ）の曲線下面積ＡＵＣに略等価になる。従って、図１０に示すように、上矢状静脈洞の時間濃度曲線Ｃｓｓｓ（ｔ）の曲線下面積ＡＵＣｓｓｓに対して、Ｓ６で作成した脳動脈の時間濃度曲線Ｃａ（ｔ）の曲線下面積ＡＵＣａが略等価になるように、Ｓ６で作成した脳動脈の時間濃度曲線Ｃａ（ｔ）の各時刻値に、ＡＵＣ（ｓｓｓ／ＡＵＣａ）を乗算することで補正する（Ｓ１０）。
【０１５３】
次に、以上のようにノイズ及びパーシャルボリューム効果が抑圧された図１３（ａ）に示す脳動脈の時間濃度曲線Ｃａ（ｔ）を使って、脳組織（毛細血管）の血流動態の様子を定量化する。そのためにまず、脳組織上の各画素について、図１３（ｂ）に示す時間濃度曲線Ｃｉ（ｔ）が作成される（Ｓ１１）。
【０１５４】
次に、Ｓ１２に示すように、左右エリアで別々の脳動脈時間濃度曲線Ｃａ（ｔ）を使って、画素ごとに、脳動脈時間濃度曲線Ｃａ（ｔ）を入力関数、脳組織の時間濃度曲線Ｃｉ（ｔ）を出力関数として、トレーサーが毛細血管を通過する過程の特徴を、伝達関数ＭＴＦとして求める。つまり、左エリアの脳組織の時間濃度曲線Ｃｉ（ｔ）に対しては同じ左エリアの脳動脈の時間濃度曲線Ｃａ（ｔ）を使い、また右エリアの脳組織の時間濃度曲線Ｃｉ（ｔ）に対しては同じ右エリアの脳動脈の時間濃度曲線Ｃａ（ｔ）を使って、伝達関数ＭＴＦを求める。また、上述したように脳動脈時間濃度曲線Ｃａ（ｔ）は、ＡＣＡ，ＭＣＡ，ＰＣＡごとに作成されるので、各Ｃａ（ｔ）ごとに伝達関数ＭＴＦの計算が繰り返される。
【０１５５】
ここでは、図１４に示すように、脳動脈の時間濃度曲線Ｃａ（ｔ）と、毛細血管の時間濃度曲線Ｃｉ（ｔ）との間で成り立つ理想的な関係を解析モデルとして用いる。つまり脳組織に入る直前の血管から造影剤（トレーサー）を瞬間的に注入した場合、毛細血管を含む脳組織単位体積（１画素）内の時間濃度曲線は立ち上がりが垂直で、若干の勾配を持って立ち下がる形になり、これが組織のインパルス応答、すなわち伝達関数を表している。そこで伝達関数を、矩形関数で近似する。これをｂｏｘ−ＭＴＦ法と呼ぶ。
【０１５６】
図１５にｂｏｘ−ＭＴＦ法の原理を示している。脳動脈の時間濃度曲線Ｃａ（ｔ）に対して、矩形関数で表される伝達関数ｂｏｘ−ＭＴＦをコンボリュートし、それにより得られたＣｉ´（ｔ）と、Ｓ１１で作成した実測Ｃｉ（ｔ）との残差を評価し、その残差の二乗和を減少させるように伝達関数ｂｏｘ−ＭＴＦを修正する。このルーチンを繰り返すことで、残差を最小化させる。
【０１５７】
残差を最小化する伝達関数ｂｏｘ−ＭＴＦに基づいて、図１４に示すように、ＣＢＰ、ＣＢＶ、ＭＴＴが計算され（Ｓ１３）、またＳ１２で最小化した残差の二乗和がそのままＥｒｒとして取り込まれる。厳密には、  
ＣＢＰ＝ＣＢＰ  
ＣＢＶ＝（１−Ｈｔ）／（１−ｂ＊Ｈｔ）＊ＣＢＶ'  
ＭＴＴ＝（１−Ｈｔ）／（１−ｂ＊Ｈｔ）＊ＭＴＴ'  
で補正される。ここで、Ｈｔは大血管のヘマトクリット値であり、ｂ＊Ｈｔは末梢血管のヘマトクリット値（一般的にはｂは０．７程度）である。
【０１５８】
こうしてｂｏｘ−ＭＴＦ法により得られた伝達関数から計算されたＣＢＰ、ＣＢＶ、ＭＴＴ、Ｅｒｒに対して、図１６に示すように、個々に出力範囲（適正範囲）が設定される（Ｓ１４）。それぞれ対応する出力範囲内の値を持つ画素については、その値が維持され、出力範囲外の値を持つ画素にはその値から、例えば表示上で黒に対応する値に置換される（Ｓ１５）。
【０１５９】
次に、図１７に示すように、出力適正化を受けたＣＢＰ、ＣＢＶ、ＭＴＴ、Ｅｒｒからそれぞれのマップが生成される（Ｓ１６）。ＣＢＰ、ＣＢＶ、ＭＴＴ、Ｅｒｒのインデックスは、前大脳動脈ＡＣＡ，中大脳動脈ＭＣＡ，後大脳動脈ＰＣＡごと、しかも各スライスに対して個々に計算される。従って、マップは、１スライスでも、図１８に示すように、４×３＝１２枚になる。設定する脳動脈数を増加すれば、その増加数の４倍でマップは増える。このような多くのマップを総合的に評価することは現実的ではない。そこで、マップ枚数を減少させるために、マップを合成する（Ｓ１７）。
【０１６０】
合成方法としては、前大脳動脈ＡＣＡのＣＢＰマップと、中大脳動脈ＭＣＡのＣＢＰマップと、後大脳動脈ＰＣＡのＣＢＰマップとを、前大脳動脈ＡＣＡの残差Ｅｒｒ、中大脳動脈ＭＣＡの残差Ｅｒｒ、後大脳動脈ＰＣＡの残差Ｅｒｒに基づいて合成する。例えば前大脳動脈ＡＣＡの時間濃度曲線Ｃａ（ｔ）とその支配下にある脳組織の時間濃度曲線Ｃｉ（ｔ）とから伝達関数ＭＴＦを求めた場合、その残差Ｅｒｒは比較的少なく、逆に、支配下にない脳組織の時間濃度曲線Ｃｉ（ｔ）から伝達関数ＭＴＦを求めた場合、その残差Ｅｒｒは比較的多くなる。つまり、残差Ｅｒｒは各脳動脈の支配可能性を表している。
【０１６１】
従って、各画素ごとに、前大脳動脈ＡＣＡのＣＢＰ値と、中大脳動脈ＭＣＡのＣＢＰ値と、後大脳動脈ＰＣＡのＣＢＰ値との中から、最も残差Ｅｒｒの低い値に対応するＣＢＰ値をその画素の値として選択する。こうして合成されたマップは、前大脳動脈ＡＣＡ，中大脳動脈ＭＣＡ，後大脳動脈ＰＣＡの支配下にある可能性の高い脳組織のＣＢＰ値から構成される。他のインデックスＣＢＶ、ＭＴＴのマップ合成についても同様である。
【０１６２】
また、最も残差Ｅｒｒが低い値を示す脳動脈が、前大脳動脈ＡＣＡと中大脳動脈ＭＣＡと後大脳動脈ＰＣＡとのいずれであるかを選択し、選択された前大脳動脈ＡＣＡ、中大脳動脈ＭＣＡ又は後大脳動脈ＰＣＡに応じて相違する色情報を与えることにより、前大脳動脈ＡＣＡの支配領域と、中大脳動脈ＭＣＡの支配領域と、後大脳動脈ＰＣＡの支配領域とを色分けにより区別された支配マップを作成するようにしても良い。
【０１６３】
ここで、マップ合成について以下に詳細に説明する。動脈に対応する位置にある画素から得られる時間濃度曲線は動脈血中造影剤濃度を反映しており、これに上記したコヒーレントレグレッション法等を適用して正確な動脈血中造影剤濃度の時間濃度曲線を得ることができる。このような脳動脈の時間濃度曲線は動脈ごとに作成可能であり、それぞれ血行状態によって違いがある。特に脳血管障害等を起こしている症例に於いてはこの違いが著しい場合がある。例えばＫ箇所の動脈で得た脳動脈の時間濃度曲線をＡｋ（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）を観察することによって、解剖学的知識なしに、どの組織がどの動脈の支配を受けているかを観察することができる。
【０１７６】
ここで、図６に戻る。図２０（ａ）、図２０（ｂ）、図２０（ｃ）、図２０（ｄ）に示すように、以上のように合成された、又は各脳動脈で単独のＣＢＰマップ、ＣＢＶマップ、ＭＴＴマップ、Ｅｒｒマップに対して、複数画素を含む関心領域ＲＯＩを設定し（Ｓ１８）、そのＲＯＩ内の画素値（ＣＢＰ値，ＣＢＶ値、ＭＴＴ値、Ｅｒｒ値）の平均値（ＣＢＰ平均値、ＣＢＶ平均値、ＭＴＴ平均値、Ｅｒｒ平均値）を計算し（Ｓ１９）、その平均値を診断材料とすることがある。この平均化に際して、上記ステップＳ１４でＣＢＰ、ＣＢＶ、ＭＴＴ、Ｅｒｒ各々に対して適正範囲を設定し、その範囲内の値を維持し、その範囲から外れた値は、例えば黒色表現に対応した最小値に置換したので、この置換した値を含めて平均化すると、その平均値には当然にして誤差が含まれてしまう。そのためこの平均化処理にあたっては、適正範囲内の値だけを選択して、または置換した値を除外して、平均化処理をすることが必要である。
【０１７７】
また、この平均化のための関心領域ＲＯＩの設定にあたっては、ＣＢＰマップ、ＣＢＶマップ、ＭＴＴマップ、Ｅｒｒマップのいずれかのマップ上で当該関心領域ＲＯＩを設定すれば、そのＲＯＩが他のマップにも共通で用いられるようになっており、、ＲＯＩ設定作業の簡素化を図り、また、同じＲＯＩに関する平均値（ＣＢＰ平均値、ＣＢＶ平均値、ＭＴＴ平均値、Ｅｒｒ平均値）の計算を可能としている。
【０１７８】
以上のように本実施形態によれば、コヒーレントフィルタ又はコヒーレントレグレッションを併用することにより、空間及び時間分解能の低下を抑えて、ノイズを抑制し、それによりＣＢＰスタディの解析精度を向上することができる。
【０１７９】
なお、本発明は上記実施形態そのままに限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で構成要素を変形して具体化できる。また、上記実施形態に開示されている複数の構成要素の適宜な組み合わせにより、種々の発明を形成できる。例えば、実施形態に示される全構成要素から幾つかの構成要素を削除してもよい。さらに、異なる実施形態にわたる構成要素を適宜組み合わせてもよい。
【図面の簡単な説明】
【０１８０】
【図１】ＣＢＰスタディの原理説明図。
【図２】本発明の実施形態に係る脳組織内毛細血管の血流動態に関するインデックス演算装置の構成を示すブロック図。
【図３】本実施形態のコヒーレントフィルタによる画像処理の説明図。
【図４】本実施形態におけるコヒーレントフィルタによるノイズ抑制処理の流れを示すフローチャート。
【図５】本実施形態におけるインデックス演算処理全体の前半部のフローチャート。
【図６】本実施形態におけるインデックス演算処理全体の後半部のフローチャート。
【図７】図５のステップＳ３の分割線の一例を示す中間調画像。
【図８】図５のステップＳ４のＡＴ，ＰＴ，ＴＴを示す図。
【図９】図５のステップＳ４のＡＴマップ，ＰＴマップ，ＴＴマップを示す図。
【図１０】図５のステップＳ１０の脳動脈の時間濃度曲線の補正に関する補足図。
【図１１】図５のステップＳ７で設定される上矢状静脈洞ＲＯＩを示す中間調画像。
【図１２】図５のステップＳ６でスライス間で共通される脳動脈ＲＯＩを示す中間調画像。
【図１３】図５のステップＳ１０、Ｓ１１で作成された脳動脈の時間濃度曲線Ｃａ（ｔ）と脳組織の時間濃度曲線Ｃｉ（ｔ）との一例を示す図。
【図１４】図６のステップＳ１２のｂｏｘ−ＭＴＦ法の原理説明図。
【図１５】図６のステップＳ１２のｂｏｘ−ＭＴＦ処理の説明図。
【図１６】図６のステップＳ１４の各インデックスの出力範囲設定画面の一例を示す図。
【図１７】図６のステップＳ１６で作成されたＣＢＰ，ＣＢＶ，ＭＴＴ，Ｅｒｒの各マップの一例を示す中間調画像。
【図１８】図６のステップＳ１６で脳動脈ごとに作成されたＣＢＰマップ，ＣＢＶマップ，ＭＴＴマップ，Ｅｒｒマップを一覧で表示した中間調画像。
【図１９】図６のステップＳ１７のマップ合成法を説明するための中間調画像。
【図２０】図６のステップＳ１９で計算された平均値の表示例を示す中間調画像。
【符号の説明】
【０１８１】
１０…ガントリ部、
２０…コンピュータ装置、
１０１…Ｘ線管、
１０１ａ…高電圧発生装置、
１０２…Ｘ線検出器、
１０３…データ収集部、
３０…画像処理装置、
１０７…画像表示部、
１０９…入力部、
１０８…制御部、
１０４…前処理部、
１０５…メモリ部、
１０６…画像再構成部、
１１０…コヒーレントフィルタ処理部、
１２０…ＣＢＰスタディ処理部、
１１１…分散値推定部、
１１２…重み関数演算部、
１１３…画素値演算部（コヒーレントフィルタ部）、
１２１…ＣＢＰスタディ処理部、
１２２…ＲＯＩ設定支援部、
１２２…時間濃度曲線作成部、
１２３…脳動脈時間濃度曲線補正部、
１２４…ＭＴＦ処理部、
１２５…インデックス計算部、
１２６…マップ作成部、
１２７…マップ合成部。
【図１】
2008100121.tif 000016
【図２】
2008100121.tif 000017
【図３】
2008100121.tif 000018
【図４】
2008100121.tif 000019
【図５】
2008100121.tif 000020
【図６】
2008100121.tif 000021
【図７】
2008100121.tif 000022
【図８】
2008100121.tif 000023
【図９】
2008100121.tif 000024
【図１０】
2008100121.tif 000025
【図１１】
2008100121.tif 000026
【図１２】
2008100121.tif 000027
【図１３】
2008100121.tif 000028
【図１４】
2008100121.tif 000029
【図１５】
2008100121.tif 000030
【図１６】
2008100121.tif 000031
【図１７】
2008100121.tif 000032
【図１８】
2008100121.tif 000033
【図１９】
2008100121.tif 000034
【図２０】
2008100121.tif 000035