JP2011050410A - 血流動態解析装置、磁気共鳴イメージング装置、およびプログラム - Google Patents

Abstract

【課題】推定誤差の小さい動脈入力関数を求めることが可能な血流動態解析装置、磁気共鳴イメージング装置、血流動態解析方法、およびプログラムを提供する。
【解決手段】ＴＴＰマップＭｔ１〜ＭｔｎおよびΔＲ_２ ^＊ピーク値マップＭｖ１〜Ｍｖｎに加えて、（ΔＲ_２ ^＊）^２データ系列ＳＤ３のピーク到達時刻ＴＴＰに関する対称性を表す対称性指標ＳのマップＭｓ１〜Ｍｓｎを作成する。その後、位置座標が共通するピクセルＰα、Ｐβ、およびＰγを検出し、検出されたピクセルＰα、Ｐβ、およびＰγの位置座標に対応付けられているΔＲ_２ ^＊データ系列ＳＤ２を加算平均することにより、動脈入力関数を算出する。
【選択図】図１６

Description

本発明は、被検体の血流の動態を解析する血流動態解析装置、磁気共鳴イメージング装置、およびプログラムに関する。

被検体の脳の血流の動態を解析するために使用されるパラメータとして、脳血流量（Cerebral Blood Flow）ＣＢＦや、脳血液量（Cerebral Blood Volume）ＣＢＶなどがある。脳血流量ＣＢＦや脳血液量ＣＢＶを算出するには、動脈の位置における磁気共鳴信号の信号強度の時間変化に基づいて、動脈入力関数を予め算出しておく必要がある。動脈では、造影剤の流入に伴って磁気共鳴信号の信号強度が急激に変化するので、信号強度の時間変化を調べることによって、動脈の位置における磁気共鳴信号の信号強度の時間変化を求めることができ、動脈入力関数を算出することが可能となる。信号強度の時間変化に基づいて動脈入力関数を算出する種々の方法が提案されている（例えば、特許文献１参照）。

特開2004-057812号公報

しかし、特許文献１の方法では、体動などが原因で、動脈とは別の組織において磁気共鳴信号の信号強度が急激に変化した場合でも、動脈の位置における信号強度の時間変化と認識してしまうことがある。したがって、動脈入力関数の推定誤差が大きくなるという問題がある。

本発明は、上記の事情に鑑み、推定誤差の小さい動脈入力関数を求めることが可能な血流動態解析装置、磁気共鳴イメージング装置、およびプログラムを提供することを目的とする。

上記の問題を解決する本発明の第１の血流動態解析装置は、
造影剤が注入された被検体の各領域から各時刻における磁気共鳴信号を収集し、収集された上記磁気共鳴信号の信号強度に基づいて上記被検体の血流の動態を解析する血流動態解析装置であって、
収集された上記磁気共鳴信号に基づいて、上記被検体の各領域ごとに、上記各時刻における磁気共鳴信号の信号強度を算出し、上記信号強度の時間変化を表す第１のデータ系列を作成する第１のデータ系列作成手段と、
上記第１のデータ系列に基づいて、上記各時刻におけるスピンの緩和速度と、上記造影剤が上記被検体の各領域に到達する前におけるスピンの緩和速度との差を算出し、上記スピンの緩和速度の差の時間変化を表す第２のデータ系列を作成する第２のデータ系列作成手段と、
上記第２のデータ系列に基づいて、上記各時刻ごとに、上記スピンの緩和速度の差に関連する値を算出し、上記スピンの緩和速度の差に関連する値の時間変化を表す第３のデータ系列を作成する第３のデータ系列作成手段と、
時間軸方向に関する上記第３のデータ系列の対称性を表す対称性指標を算出する対称性指標算出手段と、
上記対称性指標に基づいて、上記第２のデータ系列作成手段が作成したデータ系列の中から、動脈入力関数を決定するために使用される上記第２のデータ系列を選択するデータ系列選択手段と、
を有している。
上記の問題を解決する本発明の第２の血流動態解析装置は、
造影剤が注入された被検体の各領域から各時刻における磁気共鳴信号を収集し、収集された上記磁気共鳴信号の信号強度に基づいて上記被検体の血流の動態を解析する血流動態解析装置であって、
収集された上記磁気共鳴信号に基づいて、上記被検体の各領域ごとに、上記各時刻における磁気共鳴信号の信号強度を算出し、上記信号強度の時間変化を表す第１のデータ系列を作成する第１のデータ系列作成手段と、
上記第１のデータ系列に基づいて、上記各時刻における上記造影剤の濃度を算出し、上記造影剤の濃度の時間変化を表す第２のデータ系列を作成する第２のデータ系列作成手段と、
上記第２のデータ系列に基づいて、上記各時刻ごとに、上記造影剤の濃度に関連する値を算出し、上記造影剤の濃度に関連する値の時間変化を表す第３のデータ系列を作成する第３のデータ系列作成手段と、
時間軸方向に関する上記第３のデータ系列の対称性を表す対称性指標を算出する対称性指標算出手段と、
上記対称性指標に基づいて、上記第２のデータ系列作成手段が作成したデータ系列の中から、動脈入力関数を決定するために使用される上記第２のデータ系列を選択するデータ系列選択手段と、
を有している。
上記の問題を解決する本発明の磁気共鳴イメージング装置は、本発明の第１の血流動態解析装置又は本発明の第２の血流動態解析装置を有している。
上記の問題を解決する本発明の第１のプログラムは、
造影剤が注入された被検体の各領域から各時刻における磁気共鳴信号を収集し、収集された上記磁気共鳴信号の信号強度に基づいて上記被検体の血流の動態を解析するためのプログラムであって、
収集された上記磁気共鳴信号に基づいて、上記被検体の各領域ごとに、上記各時刻における磁気共鳴信号の信号強度を算出し、上記信号強度の時間変化を表す第１のデータ系列を作成する第１のデータ系列作成処理と、
上記第１のデータ系列に基づいて、上記各時刻におけるスピンの緩和速度と、上記造影剤が上記被検体の各領域に到達する前におけるスピンの緩和速度との差を算出し、上記スピンの緩和速度の差の時間変化を表す第２のデータ系列を作成する第２のデータ系列作成処理と、
上記第２のデータ系列に基づいて、上記各時刻ごとに、上記スピンの緩和速度の差に関連する値を算出し、上記スピンの緩和速度の差に関連する値の時間変化を表す第３のデータ系列を作成する第３のデータ系列作成処理と、
時間軸方向に関する上記第３のデータ系列の対称性を表す対称性指標を算出する対称性指標算出処理と、
上記対称性指標に基づいて、上記第２のデータ系列作成処理により作成されたデータ系列の中から、動脈入力関数を決定するために使用される上記第２のデータ系列を選択するデータ系列選択処理と、
を計算機に実行させるためのプログラムである。
上記の問題を解決する本発明の第２のプログラムは、
造影剤が注入された被検体の各領域から各時刻における磁気共鳴信号を収集し、収集された上記磁気共鳴信号の信号強度に基づいて上記被検体の血流の動態を解析するためのプログラムであって、
収集された上記磁気共鳴信号に基づいて、上記被検体の各領域ごとに、上記各時刻における磁気共鳴信号の信号強度を算出し、上記信号強度の時間変化を表す第１のデータ系列を作成する第１のデータ系列作成処理と、
上記第１のデータ系列に基づいて、上記各時刻における上記造影剤の濃度を算出し、上記造影剤の濃度の時間変化を表す第２のデータ系列を作成する第２のデータ系列作成処理と、
上記第２のデータ系列に基づいて、上記各時刻ごとに、上記造影剤の濃度に関連する値を算出し、上記造影剤の濃度に関連する値の時間変化を表す第３のデータ系列を作成する第３のデータ系列作成処理と、
時間軸方向に関する上記第３のデータ系列の対称性を表す対称性指標を算出する対称性指標算出処理と、
上記対称性指標に基づいて、上記第２のデータ系列作成処理により作成されたデータ系列の中から、動脈入力関数を決定するために使用される上記第２のデータ系列を選択するデータ系列選択処理と、
を計算機に実行させるためのプログラムである。

本発明では、データ系列の対称性を表す対称性指標を算出し、対称性指標に基づいて、動脈入力関数を決定するために使用されるデータ系列を選択している。したがって、推定誤差の小さい動脈入力関数を得ることができる。

本発明の第１の実施形態の磁気共鳴イメージング装置の概略図である。 ＭＲＩ装置１の処理フローを示す図である。 被検体９に設定されたスライスＡ１〜Ａｎの一例である。 スライスＡ１〜Ａｎから得られるフレーム画像を示す概念図である。 スライスＡ１〜Ａｎと、フレーム画像との対応関係を示す図である。 スライスＡｋの領域Ｒαにおける磁気共鳴信号の信号強度Ｓ（ｔ）の時間変化を表す信号強度データ系列の一例を示す図である。 スライスＡ１〜Ａｎの各領域と、作成された信号強度データ系列ＳＤ１との対応関係を示す概念図である。 スライスＡｋの領域ＲαにおけるΔＲ_２ ^＊データ系列ＳＤ２を示す図である。 スライスＡ１〜Ａｎの各領域と、作成されたΔＲ_２ ^＊データ系列ＳＤ２との対応関係を示す概念図である。 ΔＲ_２ ^＊（ｔ）のピーク値Ｖのマップ、およびピーク到達時刻ＴＴＰのマップを概略的に示す図である。 ＴＴＰマップＭｔ１〜Ｍｔｎの各ピクセルにおけるピーク到達時刻ＴＴＰと、ピクセル数との関係を表すグラフである。 ＴＴＰマップＭｔ１〜Ｍｔｎにおいて、ピーク到達時刻ＴＴＰの値が範囲Ｗ（０＜ＴＴＰ＜Ｑ１−σ）に含まれるピクセルを示す概略図である。 ΔＲ_２ ^＊ピーク値マップＭｖ１〜Ｍｖｎの各ピクセルにおけるΔＲ_２ ^＊のピーク値Ｖと、ピクセル数との関係を表すグラフである。 ΔＲ_２ ^＊ピーク値マップＭｖ１〜Ｍｖｎにおいて、ピーク値Ｖの値が範囲Ｗ（Ｖ＞Ｑ２＋σ）に含まれるピクセルを示す概略図である。 スライスＡｋの領域Ｒαにおける（ΔＲ_２ ^＊）^２データ系列ＳＤ３を示すグラフである。 対称性指標Ｓの算出方法の説明図である。 対称性指標Ｓのマップを概略的に示す図である。 対称性指標マップＭｓ１〜Ｍｓｎにおいて、対称性指標ＳがＳ＝０．９以上になるピクセルを示す概略図である。 検出されたピクセルを示す図である。 位置座標（ｘ’α、ｙ’α、ｚ’α）、（ｘ’β、ｙ’β、ｚ’β）、および（ｘγ’、ｙγ’、ｚγ’）に対応付けられているΔＲ_２ ^＊データ系列ＳＤ２を示す概念図である。 スライスＡｋの領域ＲαにおけるΔＲ_２ ^＊データ系列ＳＤ２の対称性指標Ｓ’を算出する場合の説明図である。 第２の実施形態のＭＲＩ装置の処理フローを示す図である。 作成されたＣｍデータ系列ＳＤ４を示す概念図である。 Ｃｍデータ系列ＳＤ４のピーク値Ｃmaxのマップ、およびピーク到達時刻ＴＴＰのマップを概略的に示す図である。 スライスＡｋの領域Ｒαにおける（Ｃｍ）^２データ系列ＳＤ５を示すグラフである。 対称性指標Ｓのマップを概略的に示す図である。

以下、図面を参照しながら、発明を実施するための形態を詳細に説明するが、本発明は、以下の実施形態に限定されることはない。

（１）第１の実施形態
図１は、本発明の第１の実施形態の磁気共鳴イメージング装置の概略図である。

磁気共鳴イメージング装置（以下、ＭＲＩ（Magnetic Resonance Imaging）装置と呼ぶ）１は、コイルアセンブリ２、テーブル３、造影剤注入装置４、受信コイル５と、制御装置６と、入力装置７と、表示装置８とを有している。

コイルアセンブリ２は、被検体９が収容されるボア２１と、超伝導コイル２２と、勾配コイル２３と、送信コイル２４とを有している。超伝導コイル２２は静磁場Ｂ0を印加し、勾配コイル２３は勾配パルスを印加し、送信コイル２４はＲＦパルスを送信する。

テーブル３は、被検体９を搬送するためのクレードル３１を有している。クレードル３１によって、被検体９はボア２１に搬送される。

造影剤注入装置４は、被検体９に造影剤を注入する。

受信コイル５は、被検体９の頭部９ａに取り付けられている。受信コイル５が受信したＭＲ（Magnetic Resonance）信号は、制御装置６に伝送される。

制御装置６は、シーケンサ６１〜中央処理装置６５を有している。

シーケンサ６１は、中央処理装置６５の制御を受けて、ＲＦパルスの情報（中心周波数、バンド幅など）を送信器６２に送り、勾配磁場の情報（勾配磁場の強度など）を勾配磁場電源６３に送る。

送信器６２は、シーケンサ６１から送られた情報に基づいて、ＲＦコイル２４を駆動する駆動信号を出力する。

勾配磁場電源６３は、シーケンサ６１から送られた情報に基づいて、勾配コイル２３を駆動する駆動信号を出力する。

受信器６４は、受信コイル５で受信された磁気共鳴信号を信号処理し、中央処理装置６５に伝送する。

中央処理装置６５は、シーケンサ６１および表示装置８に必要な情報を伝送したり、受信器６４から受け取った信号に基づいて画像を再構成するなど、ＭＲＩ装置１の各種の動作を実現するように、ＭＲＩ装置１の各部の動作を総括する。中央処理装置６５は、例えばコンピュータ（computer）によって構成される。尚、中央処理装置６５は、本発明における第１のデータ系列作成手段、第２のデータ系列作成手段、第３のデータ系列作成手段、対称性指標算出手段、データ系列選択手段、および特徴量算出手段の一例であり、所定のプログラムを実行することにより、これらの手段として機能する。

入力装置７は、オペレータ１０の操作に応じて、種々の命令を制御装置６に入力する。

表示装置８は種々の情報を表示する。

ＭＲＩ装置１は上記のように構成されている。次に、ＭＲＩ装置１の動作について、図２を参照しながら説明する。

図２は、ＭＲＩ装置１の処理フローを示す図である。

ステップＳ１では、被検体９の頭部９ａの造影撮影が行われる。オペレータ１０は、造影撮影を行うために、入力装置７（図１参照）を操作して、被検体９にスライスを設定する。

図３は、被検体９に設定されたスライスＡ１〜Ａｎの一例である。

被検体９には、ｎ枚のスライスＡｊ（ｊ＝１〜ｎの整数）が設定されている。スライスの枚数は、例えば、ｎ＝１２である。スライスの枚数は、必要に応じて、任意の枚数を設定することができる。

オペレータ１０は、スライスＡ１〜Ａｎを設定した後、入力装置７を操作して、被検体９を撮影する撮影命令を入力する。撮影命令が入力されると、造影剤注入装置４は被検体９に造影剤を注入するとともに、シーケンサ６１は、中央処理装置６５の制御を受けて、パルスシーケンスを実行するための情報（ＲＦパルスの情報、勾配磁場の情報）を送信器６２および勾配磁場電源６３に送る。これによって、パルスシーケンスが実行され、被検体９から磁気共鳴信号が収集される。第１の実施形態では、マルチスライススキャンにより、各スライスＡ１〜Ａｎから磁気共鳴信号が収集される。磁気共鳴信号は、受信コイル５で受信され、受信器６４で処理され、中央処理装置６５に伝送される。

被検体９から磁気共鳴信号を収集した後、ステップＳ２に進む。

ステップＳ２では、収集した磁気共鳴信号に基づいて、磁気共鳴信号の信号強度Ｓ（ｔ）の時間変化を表す信号強度データ系列を作成する。以下に、信号強度データ系列を作成する方法について説明する。

図４は、スライスＡ１〜Ａｎから得られるフレーム画像を示す概念図である。

第１の実施形態では、各スライスＡ１〜Ａｎからｍ枚のフレーム画像が得られる。したがって、１枚のスライスにつき、ｍ枚のフレーム画像が得られる。例えば、フレーム画像の枚数ｍ＝８５枚である。

図４は、被検体９の頭部９ａに設定されたｎ枚のスライスＡ１〜Ａｎにおけるフレーム画像を、収集順序に従って時系列に並べて示している。図４では、スライスＡ１〜ＡｎのうちのスライスＡ１、Ａｋ、およびＡｎにおけるフレーム画像について、それぞれ符号Ｆ１１〜Ｆ１ｍ、Ｆｋ１〜Ｆｋｍ、およびＦｎ１〜Ｆｎｍで示されている。例えば、スライスＡｋにおけるフレーム画像Ｆｋ１〜Ｆｋｍは、時刻ｔｋ１〜ｔｋｍにおいて収集された画像である。

図５は、スライスＡ１〜Ａｎと、フレーム画像との対応関係を示す図である。

図５（ａ）は、スライスＡ１〜Ａｎを示す図であり、図５（ｂ）は、各スライスＡ１〜Ａｎから収集されたフレーム画像を、スライス別に示す図である。

スライスＡｊ（ｊ＝１〜ｎの整数）の領域Ｒｉ（ｉ＝１〜ｚの整数）の位置座標は、ｘ’軸、ｙ’軸、およびｚ’軸によって規定されている。例えば、スライスＡｋ内の領域Ｒαの位置座標は、（ｘ’α、ｙ’α、ｚ’α）である。

スライスＡｊ（ｊ＝１〜ｎの整数）から収集されたフレーム画像の各ピクセルは、スライスＡｊの各領域Ｒｉから収集した磁気共鳴信号の信号強度Ｓ(ｔ)を表している。例えば、スライスＡｋから収集されたフレーム画像Ｆｋ１〜ＦｋｍのピクセルＰαは、スライスＡｋの領域Ｒαから収集した磁気共鳴信号の信号強度Ｓ(ｔ)を表している。

ステップＳ２では、収集した磁気共鳴信号に基づいて、スライスＡｊ（ｊ＝１〜ｎの整数）の領域Ｒｉ（ｉ＝１〜ｚの整数）ごとに、磁気共鳴信号の信号強度Ｓ（ｔ）の時間変化を表す信号強度データ系列を作成する。

図６は、スライスＡｋの領域Ｒαにおける磁気共鳴信号の信号強度Ｓ（ｔ）の時間変化を表す信号強度データ系列の一例を示す図である。

図６（ａ）は、スライスＡｋのフレーム画像Ｆｋ１〜Ｆｋｍを示しており、図６（ｂ）は、磁気共鳴信号の信号強度Ｓ（ｔ）の時間変化を表す信号強度データ系列ＳＤ１を示している。

信号強度データ系列ＳＤ１は、時系列的に並ぶｍ個のデータＤ１〜Ｄｍを有している。各データＤ１〜Ｄｍは、フレーム画像Ｆｋ１〜ＦｋｍのピクセルＰαの信号強度Ｓ（ｔ）を表している。例えば、データＤ１は、フレーム画像Ｆｋ１のピクセルＰαの信号強度Ｓ（ｔ）を表しており、データＤｍは、フレーム画像ＦｋｍのピクセルＰαの信号強度Ｓ（ｔ）を表している。信号強度データ系列ＳＤ１を参照すると、スライスＡｋの領域Ｒαにおいては、信号強度Ｓ（ｔ）は、時刻ｔｓを境にして急激に減少するが、時間の経過かと共に再び増大している。領域Ｒαにおける信号強度Ｓ（ｔ）は、造影剤が領域Ｒαに到達すると急激に小さい値になる。したがって、時刻ｔｓは、領域Ｒαに造影剤が到達した時刻と考えられる。

また、信号強度データ系列ＳＤ１から、信号強度Ｓ（ｔ）のベースラインＳ０も算出する。ベースラインＳ０は、造影剤が領域Ｒαに到達する前における信号強度を表している。ベースラインＳ０は、例えば、時刻ｔｓよりも前の信号強度Ｓ（ｔ）の平均値として算出することができる。ただし、時刻ｔｋ１における信号強度Ｓ（ｔ）の値は、他の時刻における信号強度Ｓ（ｔ）よりも極端に大きいので、ベースラインＳ０を算出するために用いる信号強度Ｓ（ｔ）としては適していないと考えられる。したがって、時刻ｔｋ１における信号強度Ｓ（ｔ）の値は、ベースラインＳ０を算出するために用いる信号強度Ｓ（ｔ）から除外する。

図６には、スライスＡｋの領域Ｒαにおける信号強度データ系列ＳＤ１が示されているが、スライスＡｋの各領域Ｒｊごとに、信号強度Ｓ（ｔ）の時間変化を表す信号強度データ系列ＳＤ１を作成する。

尚、図６では、スライスＡｋの信号強度データ系列ＳＤ１について説明されている。しかし、スライスＡｋ以外の他のスライスＡｊの各領域Ｒｊについても、信号強度Ｓ（ｔ）の時間変化を表す信号強度データ系列ＳＤ１を作成する。したがって、スライスＡ１〜Ａｎの各領域ごとに、信号強度データ系列ＳＤ１が作成される。作成された信号強度データ系列ＳＤ１は、スライスＡ１〜Ａｎの各領域に対応付けられる（図７参照）。

図７は、スライスＡ１〜Ａｎの各領域と、作成された信号強度データ系列ＳＤ１との対応関係を示す概念図である。

スライスＡ１〜Ａｎの各領域には、作成された信号強度データ系列ＳＤ１が対応付けられている。例えば、位置座標（ｘ’α、ｙ’α、ｚ’α）における領域Ｒαには、図６（ｂ）に示す信号強度データ系列ＳＤ１が対応付けられている。したがって、位置座標が決まれば、対応する信号強度データ系列ＳＤ１を取り出すことができる。

信号強度データ系列ＳＤ１を作成した後、ステップＳ３に進む。

ステップＳ３では、信号強度データ系列ＳＤ１に基づいて、スライスＡｊ（ｊ＝１〜ｎ）の領域Ｒｉ（ｉ＝１〜ｚ）ごとに、ΔＲ_２ ^＊（ｔ）の時間変化を表すΔＲ_２ ^＊データ系列ＳＤ２を作成する。ΔＲ_２ ^＊（ｔ）は、以下の式で表される。
ΔＲ_２ ^＊（ｔ）＝Ｒ_２ ^＊（ｔ）−Ｒ_２ ^＊（ｔ_０）・・・（１）
ここで、Ｒ_２ ^＊（ｔ）：スライスＡｊの領域Ｒｉにおける時刻ｔのスピンの緩和速度
Ｒ_２ ^＊（ｔ_０））：造影剤がスライスＡｊの領域Ｒｉに到達する前のスピンの緩和速度

したがって、式（１）から、ΔＲ_２ ^＊（ｔ）は、時刻ｔのスピンの緩和速度と、造影剤が到達する前のスピンの緩和速度との差を表していることがわかる。式（１）の右辺の項Ｒ_２ ^＊（ｔ）およびＲ_２ ^＊（ｔ_０）は、以下の式に従って算出することができる。
Ｒ_２ ^＊（ｔ） ＝−（１／ＴＥ）ｌｏｇ（Ｓ（ｔ）） ・・・（２）
Ｒ_２ ^＊（ｔ_０）＝−（１／ＴＥ）ｌｏｇ（Ｓ０） ・・・（３）
ここで、ＴＥ：エコー時間、Ｓ（ｔ）：時刻ｔにおける信号強度、Ｓ０：信号強度のベースライン

式（２）および（３）を式（１）に代入し、整理すると、以下の式が得られる。
ΔＲ_２ ^＊（ｔ）＝−（１／ＴＥ）ｌｏｇ（Ｓ（ｔ）／Ｓ０）・・・（４）

エコー時間ＴＥは、予め求められる値であり、信号強度Ｓ（ｔ）およびベースラインＳ０は、ステップＳ２において算出された値である。したがって、エコー時間ＴＥ、信号強度Ｓ（ｔ）、およびベースラインＳ０は、全て既知の値であるので、ＴＥ、Ｓ（ｔ）、およびＳ０を式（４）に代入することによって、緩和速度の差ΔＲ_２ ^＊（ｔ）の時間変化を表すΔＲ_２ ^＊データ系列ＳＤ２を作成することができる。

図８は、スライスＡｋの領域ＲαにおけるΔＲ_２ ^＊データ系列ＳＤ２を示す図である。

図８には、スライスＡｋの領域Ｒαにおける信号強度データ系列ＳＤ１も示されている。

ΔＲ_２ ^＊データ系列ＳＤ２を参照すると、緩和速度の差ΔＲ_２ ^＊（ｔ）は、時刻ｔｓを境にして急激に大きくなり、時刻ｔｐにおいて最大値となるが、再び減少することが分かる。

尚、図８では、スライスＡｋの領域ＲαにおけるΔＲ_２ ^＊データ系列ＳＤ２が示されている。しかし、スライスＡｋの他の領域についてもΔＲ_２ ^＊データ系列ＳＤ２を作成し、スライスＡｋ以外の他のスライスＡｊの各領域についても、ΔＲ_２ ^＊データ系列ＳＤ２を作成する。したがって、スライスＡ１〜Ａｎの各領域ごとに、ΔＲ_２ ^＊データ系列ＳＤ２が作成される。作成されたΔＲ_２ ^＊データ系列ＳＤ２は、スライスＡ１〜Ａｎの各領域ごとに対応付けられる（図９参照）。

図９は、スライスＡ１〜Ａｎの各領域と、作成されたΔＲ_２ ^＊データ系列ＳＤ２との対応関係を示す概念図である。

スライスＡ１〜Ａｎの各領域には、作成されたΔＲ_２ ^＊データ系列ＳＤ２が対応付けられている。例えば、位置座標（ｘ’α、ｙ’α、ｚ’α）における領域Ｒαには、図８に示すΔＲ_２ ^＊データ系列ＳＤ２が対応付けられている。したがって、位置座標が決まれば、対応するΔＲ_２ ^＊データ系列ＳＤ２を取り出すことができる。尚、ステップＳ２において既に信号強度データ系列ＳＤ１が対応付けられているので、位置座標ごとに、信号強度データ系列ＳＤ１とΔＲ_２ ^＊データ系列ＳＤ２とが対応付けられている。

ΔＲ_２ ^＊データ系列ＳＤ２を作成した後、ステップＳ４に進む。

ステップＳ４では、ΔＲ_２ ^＊データ系列ＳＤ２の特徴量を算出し、特徴量のマップを作成する。第１の実施形態では、ΔＲ_２ ^＊データ系列ＳＤ２の特徴量として、以下の２個の特徴量を算出する。
（１）ΔＲ_２ ^＊データ系列ＳＤ２のピーク部ＰＥにおけるピーク値Ｖ
（２）ピーク値Ｖにおける時刻（ピーク到達時刻）ＴＴＰ

図９には、スライスＡｋの領域ＲαにおけるΔＲ_２ ^＊データ系列ＳＤ２のピーク値Ｖ（＝Ｖα）およびピーク到達時刻ＴＴＰ（＝ＴＴＰα）が示されている。しかし、スライスＡｋの他の領域におけるΔＲ_２ ^＊データ系列ＳＤ２のピーク値Ｖおよびピーク到達時刻ＴＴＰも算出される。更に、スライスＡｋ以外の他のスライスの各領域におけるΔＲ_２ ^＊データ系列ＳＤ２のピーク値Ｖおよびピーク到達時刻ＴＴＰも算出される。

尚、図９では、ΔＲ_２ ^＊データ系列ＳＤ２の中で、ΔＲ_２ ^＊（ｔ）の最大値をピーク値Ｖαとしている。しかし、ΔＲ_２ ^＊データ系列ＳＤ２に最もよく当てはまるフィッティング曲線を求め、フィッティング曲線の最大値をピーク値Ｖαとするなど、別の方法でピーク値Ｖαを算出してもよい。

図１０は、ΔＲ_２ ^＊データ系列ＳＤ２のピーク値Ｖのマップ、およびΔＲ_２ ^＊データ系列ＳＤ２のピーク到達時刻ＴＴＰのマップを概略的に示す図である。

ピーク到達時刻ＴＴＰのマップ（以下、「ＴＴＰマップ」と呼ぶ）Ｍｔ１〜Ｍｔｎの各ピクセルのピクセル値は、スライスＡ１〜Ａｎの各領域におけるピーク到達時刻ＴＴＰの値を表している。例えば、ＴＴＰマップＭｔｋのピクセルＰαのピクセル値は、スライスＡｋの領域Ｒαにおけるピーク到達時刻ＴＴＰα（図９参照）を表している。

また、ピーク値Ｖのマップ（以下、「ΔＲ_２ ^＊ピーク値マップ」と呼ぶ）Ｍｖ１〜Ｍｖｎの各ピクセルのピクセル値は、スライスＡ１〜Ａｎの各領域におけるピーク値Ｖの値を表している。例えば、ΔＲ_２ ^＊ピーク値マップＭｖｋのピクセルＰαのピクセル値は、スライスＡｋの領域Ｒαにおけるピーク値Ｖα（図９参照）を表している。

ＴＴＰマップＭｔ１〜ＭｔｎおよびΔＲ_２ ^＊ピーク値マップＭｖ１〜Ｍｖｎを作成した後、ステップＳ５に進む。

ステップＳ５では、ＴＴＰマップＭｔ１〜Ｍｔｎのピクセルの中から、動脈のピーク到達時刻ＴＴＰを表していると考えられるピクセルを抽出する。ピクセルの抽出は、以下のようにして行う。

動脈のピーク到達時刻ＴＴＰを表していると考えられるピクセルを抽出するために、先ず、ＴＴＰマップＭｔ１〜Ｍｔｎ（図１０参照）の各ピクセルにおけるピーク到達時間ＴＴＰと、ピクセル数との関係を求める。

図１１は、ＴＴＰマップＭｔ１〜Ｍｔｎの各ピクセルにおけるピーク到達時刻ＴＴＰと、ピクセル数との関係を表すグラフである。

グラフの横軸は、ピーク到達時刻ＴＴＰの値であり、グラフの縦軸は、ＴＴＰマップＭｔ１〜Ｍｔｎに含まれるピクセルの中で、ピーク到達時刻ＴＴＰが同じ値になるピクセルの数ＮＰである。

グラフには、ピーク到達時刻ＴＴＰとピクセル数ＮＰとの関係を表す曲線Ｃｔが示されている。

動脈のピーク到達時刻ＴＴＰは、一般的に小さい値になる傾向がある。したがって、ピーク到達時刻ＴＴＰの値の小さいピクセルは、動脈のピーク到達時刻ＴＴＰを表す可能性が高いと考えられる。そこで、曲線Ｃｔに基づいて、ピーク到達時刻ＴＴＰが小さい値になる範囲Ｗを決定する。ピーク到達時刻ＴＴＰが小さい値になる範囲Ｗは、例えば、ピクセル数ＮＰが最大値ＮＰmaxになるときのピーク到達時刻ＴＴＰ（＝Ｑ１）と、ピクセル数ＮＰの標準偏差σとを算出し、０＜ＴＴＰ＜Ｑ１−σとなる範囲として決定することができる。

範囲Ｗを決定した後、ＴＴＰマップＭｔ１〜Ｍｔｎのピクセルの中から、ピーク到達時刻ＴＴＰの値が範囲Ｗ（０＜ＴＴＰ＜Ｑ１−σ）に含まれるピクセルを、動脈のピーク到達時刻ＴＴＰを表したピクセルとして抽出する（図１２参照）。

図１２は、ＴＴＰマップＭｔ１〜Ｍｔｎにおいて、ピーク到達時刻ＴＴＰの値が範囲Ｗ（０＜ＴＴＰ＜Ｑ１−σ）に含まれるピクセルを示す概略図である。

図１２では、ＴＴＰマップＭｔ１〜Ｍｔｎのピクセルの中で、ピーク到達時刻ＴＴＰの値が範囲Ｗ（０＜ＴＴＰ＜Ｑ１−σ）に含まれるピクセルを、白抜きの四角形で示してある。

ＴＴＰマップＭｔ１〜Ｍｔｎからピクセルを抽出した後、ΔＲ_２ ^＊ピーク値マップＭｖ１〜Ｍｖｎのピクセルの中から、動脈におけるΔＲ_２ ^＊のピーク値Ｖを表していると考えられるピクセルを抽出する。ピクセルの抽出は、以下のようにして行う。

動脈におけるΔＲ_２ ^＊のピーク値Ｖを表していると考えられるピクセルを抽出するために、先ず、ΔＲ_２ ^＊ピーク値マップＭｖ１〜Ｍｖｎの各ピクセルにおけるΔＲ_２ ^＊のピーク値Ｖと、ピクセル数との関係を求める。

図１３は、ΔＲ_２ ^＊ピーク値マップＭｖ１〜Ｍｖｎの各ピクセルにおけるΔＲ_２ ^＊のピーク値Ｖと、ピクセル数との関係を表すグラフである。

グラフの横軸は、ΔＲ_２ ^＊のピーク値Ｖの値であり、グラフの縦軸は、ΔＲ_２ ^＊ピーク値マップＭｖ１〜Ｍｖｎに含まれるピクセルの中で、ピーク値Ｖが同じ値になるピクセルの数ＮＰである。

グラフには、ピーク値Ｖとピクセル数ＮＰとの関係を表す曲線Ｃｖが示されている。

動脈におけるΔＲ_２ ^＊のピーク値Ｖは、一般的に大きい値になる傾向がある。したがって、ピーク値Ｖの値の大きいピクセルは、動脈のピーク値Ｖを表す可能性が高いと考えられる。そこで、曲線Ｃｖに基づいて、ピーク値Ｖが大きい値になる範囲Ｗを決定する。ピーク値Ｖが大きい値になる範囲Ｗは、例えば、ピクセル数ＮＰが最大値ＮＰmaxになるときのΔＲ_２ ^＊のピーク値Ｖ（＝Ｑ２）と、ピクセル数ＮＰの標準偏差σとを算出し、Ｖ＞Ｑ２＋σのとなる範囲として決定することができる。

範囲Ｗを決定した後、ΔＲ_２ ^＊ピーク値マップＭｖ１〜Ｍｖｎの各ピクセルの中から、ΔＲ_２ ^＊のピーク値Ｖの値が範囲Ｗ（Ｖ＞Ｑ２＋σ）に含まれるピクセルを、動脈におけるΔＲ_２ ^＊のピーク値Ｖを表したピクセルとして抽出する。

図１４は、ΔＲ_２ ^＊ピーク値マップＭｖ１〜Ｍｖｎにおいて、ピーク値Ｖの値が範囲Ｗ（Ｖ＞Ｑ２＋σ）に含まれるピクセルを示す概略図である。

図１４では、ΔＲ_２ ^＊ピーク値マップＭｖ１〜Ｍｖｎのピクセルの中で、ΔＲ_２ ^＊のピーク値Ｖが範囲Ｗ（Ｖ＞Ｑ２＋σ）に含まれるピクセルを、白抜きの四角形で示してある。

ΔＲ_２ ^＊ピーク値マップＭｖ１〜Ｍｖｎからピクセルを抽出した後、ステップＳ６に進む。

ステップＳ６では、ΔＲ_２ ^＊データ系列ＳＤ２に基づいて、ΔＲ_２ ^＊（ｔ）の２乗値｛ΔＲ_２ ^＊（ｔ）｝^２の時間変化を表す（ΔＲ_２ ^＊）^２データ系列ＳＤ３を作成する（図１５参照）。

図１５は、スライスＡｋの領域Ｒαにおける（ΔＲ_２ ^＊）^２データ系列ＳＤ３を示すグラフである。

図１５には、スライスＡｋの領域ＲαにおけるΔＲ_２ ^＊データ系列ＳＤ２（図９参照）と、ΔＲ_２ ^＊（ｔ）の２乗値｛ΔＲ_２ ^＊（ｔ）｝^２の時間変化を表す（ΔＲ_２ ^＊）^２データ系列ＳＤ３が示されている。

ΔＲ_２ ^＊データ系列ＳＤ２では、時刻によってはΔＲ_２ ^＊（ｔ）が負の値になることがある。例えば、時刻ｔｎ以降は、ΔＲ_２ ^＊（ｔ）＜０になっている。しかし、ΔＲ_２ ^＊（ｔ）を２乗することによって、｛ΔＲ_２ ^＊（ｔ）｝^２の値は、時刻に関わらずに正の値になる。

図１５では、スライスＡｋの領域Ｒαにおける（ΔＲ_２ ^＊）^２データ系列ＳＤ３が示されている。しかし、スライスＡｋの他の領域における（ΔＲ_２ ^＊）^２データ系列ＳＤ３も作成し、スライスＡｋ以外の他のスライスの各領域における（ΔＲ_２ ^＊）^２データ系列ＳＤ３も作成する。尚、（ΔＲ_２ ^＊）^２データ系列ＳＤ３を作成する理由については後述する。

（ΔＲ_２ ^＊）^２データ系列ＳＤ３を作成した後、ステップＳ７に進む。

ステップＳ７では、ピーク到達時刻ＴＴＰに関する（ΔＲ_２ ^＊）^２データ系列ＳＤ３の対称性を表す対称性指標Ｓを算出する。以下に、対称性指標Ｓの算出方法について、図１６を参照しながら説明する。尚、図１６では、スライスＡｋの領域Ｒαにおける（ΔＲ_２ ^＊）^２データ系列ＳＤ３（図１５（ｂ）参照）を取り上げて、対称性指標Ｓを算出するが、他の領域における（ΔＲ_２ ^＊）^２データ系列ＳＤ３の対称性指標Ｓも同様のやり方で算出することができる。

図１６は、対称性指標Ｓの算出方法の説明図である。

先ず、ピーク到達時刻ＴＴＰαよりも早い時刻に現れるデータ群ＤＧ１によって規定される領域の重心Ｇ１の時刻ｔg1を算出する。第１の実施形態では、重心Ｇ１の時刻ｔg1は、以下の式で算出される。

ただし、図１６において、時刻ｔｋ１における｛ΔＲ_２ ^＊（ｔ）｝^２は、次の時刻ｔｋ２における｛ΔＲ_２ ^＊（ｔ）｝^２よりも極端に大きいので、重心Ｇ１の時刻ｔg1を算出するために用いる｛ΔＲ_２ ^＊（ｔ）｝^２としては適していないと考えられる。したがって、時刻ｔｋ１における｛ΔＲ_２ ^＊（ｔ）｝^２の値は、重心Ｇ１の時刻ｔg1を算出するために用いる｛ΔＲ_２ ^＊（ｔ）｝^２から除外する。

重心Ｇ１の時刻ｔg1を算出した後、ピーク到達時刻ＴＴＰαと時刻ｔg1との時間差ａ１を算出する。
ａ１＝ＴＴＰα−ｔg1 ・・・（６）

また、ピーク到達時刻ＴＴＰαよりも遅い時刻に現れるデータ群ＤＧ２によって規定される領域の重心Ｇ２の時刻ｔg2を算出する。第１の実施形態では、重心Ｇ２の時刻ｔg2は、以下の式で算出される。

重心Ｇ２の時刻ｔg2を算出した後、時刻ｔg2とピーク到達時刻ＴＴＰαとの時間差ａ２を算出する。
ａ２＝ｔg2−ＴＴＰα・・・（８）

時間差ａ１およびａ２を算出した後、対称性指標Ｓを算出する。対称性指標Ｓは、以下の式に従って算出される。
Ｓ＝ｅｘｐ（−｜ａ１−ａ２｜／ａ１） ・・・（９）

一般的に、ａ１とａ２との差が小さければ小さいほど、データ群ＤＧ１とデータ群ＤＧ２は、ピーク到達時刻ＴＴＰαに関して対称的に分布していると考えられる。したがって、ａ１とａ２との差が小さければ小さいほど、ピーク到達時刻ＴＴＰαに関する（ΔＲ_２ ^＊）^２データ系列ＳＤ３の対称性は良好であると考えられる。式（９）より、ａ１とａ２との差が小さければ小さいほど、対称性指標Ｓは１に近づくので、対称性指標Ｓが１に近い値になることは、ピーク到達時刻ＴＴＰαに関する（ΔＲ_２ ^＊）^２データ系列ＳＤ３の対称性は良好であることを意味する。一方、ａ１とａ２との差が大きい場合は、対称性指標Ｓは０に近づくので、ピーク到達時刻ＴＴＰαに関するデータ系列ＳＤ３の対称性は悪いことを意味する。

尚、図１６では、スライスＡｋの領域Ｒαにおける（ΔＲ_２ ^＊）^２データ系列ＳＤ３の対称性指標Ｓを算出する例が示されている。しかし、スライスＡｋの他の領域における（ΔＲ_２ ^＊）^２データ系列ＳＤ３の対称性指標Ｓも算出し、スライスＡｋ以外の他のスライスの各領域における（ΔＲ_２ ^＊）^２データ系列ＳＤ３の対称性指標Ｓも算出する。

図１７は、対称性指標Ｓのマップを概略的に示す図である。

対称性指標Ｓのマップ（以下、「対称性指標マップ」と呼ぶ）Ｍｓ１〜Ｍｓｎの各ピクセルのピクセル値は、スライスＡ１〜Ａｎの各領域における対称性指標Ｓの値を表している。例えば、対称性指標マップＭｓｋのピクセルＰαのピクセル値は、スライスＡｋの領域Ｒαにおける対称性指標Ｓを表している。

対称性指標マップＭｓ１〜Ｍｓｎを作成した後、ステップＳ８に進む。

ステップＳ８では、対称性指標マップＭｓ１〜Ｍｓｎの中から、動脈における対称性指標Ｓを表していると考えられるピクセルを抽出する。動脈における対称性指標Ｓは、一般的に対称性指標Ｓが１に近くなる傾向があるので、対称性指標Ｓの値がＳ＝１に近いピクセルを抽出する。第１の実施形態では、対称性指標ＳがＳ＝０．９以上になるピクセルを、動脈の対称性指標Ｓを表していると考えられるピクセルとして抽出する。

図１８は、対称性指標マップＭｓ１〜Ｍｓｎにおいて、対称性指標ＳがＳ＝０．９以上になるピクセルを示す概略図である。

図１８では、対称性指標マップＭｓ１〜Ｍｓｎの中において、対称性指標ＳがＳ＝０．９以上となるピクセルを、白抜きの四角形で示してある。

ピクセルを抽出した後、ステップＳ９に進む。

ステップＳ９では、スライスＡ１〜Ａｎの各領域に対応付けられたΔＲ_２ ^＊データ系列ＳＤ２（図９参照）の中から、動脈入力関数を決定するために使用されるΔＲ_２ ^＊データ系列ＳＤ２を選択する。データ系列の選択は、以下のようにして行う。

先ず、ステップＳ９１において、ＴＴＰマップ、ΔＲ_２ ^＊ピーク値マップ、および対称性指標マップにおいて抽出されたピクセルの中から（図１８参照）、位置座標が共通するピクセルを検出する。

図１９は、検出されたピクセルを示す図である。

図１９では、検出されたピクセルは黒色の四角形で示されている。図１９を参照すると、スライスＡｋにおけるＴＴＰマップＭｔｋ、ΔＲ_２ ^＊ピーク値マップＭｖｋ、および対称性指標マップＭｓｋでは、ピクセルＰαが、共通する位置座標（ｘ’α、ｙ’α、ｚ’α）を有している。したがって、スライスＡｋにおけるＴＴＰマップＭｔｋ、ΔＲ_２ ^＊ピーク値マップＭｖｋ、および対称性指標マップＭｓｋに関しては、ピクセルＰαを検出する。同様に、他のスライスのマップについても、位置座標が同じピクセルを検出する。第１の実施形態では、ピクセルＰαの他に、ピクセルＰβおよびＰγが検出される。

ピクセルを検出した後、ステップＳ９２に進む。

ステップＳ９２では、検出されたピクセルＰα、Ｐβ、およびＰγの位置座標に対応付けられているΔＲ_２ ^＊データ系列ＳＤ２を取り出す。第１の実施形態では、検出されたピクセルＰα、Ｐβ、およびＰγの位置座標は、それぞれ（ｘ’α、ｙ’α、ｚ’α）、（ｘ’β、ｙ’β、ｚ’β）、および（ｘγ’、ｙγ’、ｚγ’）である（図１９参照）。したがって、位置座標（ｘ’α、ｙ’α、ｚ’α）、（ｘ’β、ｙ’β、ｚ’β）、および（ｘγ’、ｙγ’、ｚγ’）に対応付けられているΔＲ_２ ^＊データ系列ＳＤ２を取り出す（図２０参照）。

図２０は、位置座標（ｘ’α、ｙ’α、ｚ’α）、（ｘ’β、ｙ’β、ｚ’β）、および（ｘγ’、ｙγ’、ｚγ’）に対応付けられているΔＲ_２ ^＊データ系列ＳＤ２を示す概念図である。

これらのΔＲ_２ ^＊データ系列ＳＤ２を取り出した後、ステップＳ１０に進む。

ステップＳ１０では、取り出されたΔＲ_２ ^＊データ系列ＳＤ２を加算平均することにより、動脈入力関数を算出する。算出された動脈入力関数は、脳血流量ＣＢＦや脳血液量ＣＢＶなどを算出するために使用される。尚、第１の実施形態では、取り出されたΔＲ_２ ^＊データ系列ＳＤ２を加算平均することにより動脈入力関数を算出しているが、取り出されたΔＲ_２ ^＊データ系列ＳＤ２を重み付け平均するなど、別の方法で動脈入力関数を算出してもよい。また、取り出されたΔＲ_２ ^＊データ系列ＳＤ２が一つしかない場合は、取り出されたΔＲ_２ ^＊データ系列ＳＤ２が動脈入力関数となる。

以上説明したように、第１の実施形態では、（ΔＲ_２ ^＊）^２データ系列ＳＤ３の対称性指標Ｓを算出し、対称性指標Ｓの値に基づいて、動脈入力関数を算出するために使用されるΔＲ_２ ^＊データ系列ＳＤ２が決定される。動脈を含む領域では、対称性指標Ｓが大きくなる傾向があるので、対称性指標Ｓの値に基づいて、ΔＲ_２ ^＊データ系列ＳＤ２を選択することによって、動脈入力関数を決定するのに適したΔＲ_２ ^＊データ系列ＳＤ２を選択することができる。

尚、第１の実施形態では、（ΔＲ_２ ^＊）^２データ系列ＳＤ３の対称性指標Ｓに基づいて、動脈入力関数を決定するために使用されるΔＲ_２ ^＊データ系列ＳＤ２を選択している。これに対し、（ΔＲ_２ ^＊）^２データ系列ＳＤ３の対称性指標Ｓの代わりに、ΔＲ_２ ^＊データ系列ＳＤ２の対称性指標Ｓ’に基づいて、動脈入力関数を決定するために使用されるΔＲ_２ ^＊データ系列ＳＤ２を選択することも考えられる。しかし、ΔＲ_２ ^＊データ系列ＳＤ２の対称性指標Ｓ’に基づいて、ΔＲ_２ ^＊データ系列ＳＤ２を選択する場合、動脈入力関数を決定するのに適したΔＲ_２ ^＊データ系列ＳＤ２を必ずしも選択することができない恐れがある。以下にこの理由について、図２１を参照しながら説明する。

図２１は、スライスＡｋの領域ＲαにおけるΔＲ_２ ^＊データ系列ＳＤ２の対称性指標Ｓ’を算出する場合の説明図である。

ΔＲ_２ ^＊データ系列ＳＤ２に関する対称性指標Ｓ’を算出する場合、対称性指標Ｓ’は、以下の式に従って算出される。
Ｓ’＝ｅｘｐ（−｜ａ１’−ａ２’｜／ａ１’） ・・・（１０）
ただし、

式（１０）〜（１２）を用いて対称性指標Ｓ’を算出する場合、式（１１）および（１２）の積分範囲において、ΔＲ_２ ^＊（ｔ）が正の値を有する場合積分値は正になるが、ΔＲ_２ ^＊（ｔ）が負の値の場合積分値は負になる（例えば、ｔｎ＜ｔの範囲参照）。したがって、負の積分値が正の積分値を打ち消すことになるので、時間差ａ１’およびａ２’の誤差が大きくなる恐れがある。このため、対称性指標Ｓ’の信頼性が低くなるので、動脈入力関数を決定するのに適したΔＲ_２ ^＊データ系列ＳＤ２を必ずしも選択することができない恐れがある。

そこで、第１の実施形態では、（ΔＲ_２ ^＊）^２データ系列ＳＤ３の対称性指標Ｓに基づいて、動脈入力関数を決定するために使用されるΔＲ_２ ^＊データ系列ＳＤ２を選択している。｛ΔＲ_２ ^＊（ｔ）｝^２は、ΔＲ_２ ^＊（ｔ）の２乗値であるので、ΔＲ_２ ^＊（ｔ）が負の値であっても、｛ΔＲ_２ ^＊（ｔ）｝^２は必ず正の値になる（図１６参照）。したがって、負の積分値が正の積分値を打ち消すことがないので、誤差の少ない時間差ａ１およびａ２を算出することができ、動脈入力関数を決定するのに適したΔＲ_２ ^＊データ系列ＳＤ２を選択することができる。

尚、第１の実施形態では、（ΔＲ_２ ^＊）^２データ系列ＳＤ３の対称性指標Ｓに基づいて動脈入力関数を決定するためのΔＲ_２ ^＊データ系列ＳＤ２を選択している。しかし、｛ΔＲ_２ ^＊（ｔ）｝^ｐ（ｐは３以上の偶数）のデータ系列の対称性指標Ｓ’’を算出し、対称性指標Ｓ’’に基づいて動脈入力関数を決定するためのΔＲ_２ ^＊データ系列ＳＤ２を選択してもよい。｛ΔＲ_２ ^＊（ｔ）｝^ｐは、｛ΔＲ_２ ^＊（ｔ）｝^２と同様に必ず正の値になるので、｛ΔＲ_２ ^＊（ｔ）｝^ｐのデータ系列の対称性指標Ｓ’’を用いても、動脈入力関数を決定するのに適したΔＲ_２ ^＊データ系列ＳＤ２を選択することができる。

また、｛ΔＲ_２ ^＊（ｔ）｝^ｐのデータ系列の対称性指標Ｓ’’の代わりに、｜ΔＲ_２ ^＊（ｔ）^ｑ｜（ｑは１以上の奇数）のデータ系列の対称性指標Ｓ’’’を用いてもよい。｜ΔＲ_２ ^＊（ｔ）^ｑ｜はΔＲ_２ ^＊（ｔ）^ｑの絶対値であるので、必ず正の値になる。したがって、｜ΔＲ_２ ^＊（ｔ）^ｑ｜のデータ系列の対称性指標を用いても、動脈入力関数を決定するのに適したΔＲ_２ ^＊データ系列ＳＤ２を選択することができる。

尚、第１の実施形態では、ＴＴＰマップ、ΔＲ_２ ^＊ピーク値マップ、および対称性指標マップの順にピクセルを抽出している。しかし、別の抽出順序でピクセルを抽出してもよい。

また、第１の実施形態では、ステップＳ４において、ΔＲ_２ ^＊データ系列ＳＤ２の特徴量として、ピーク到達時刻ＴＴＰおよびΔＲ_２ ^＊ピーク値Ｖを算出している。しかし、別の特徴量を算出してもよい。

また、第１の実施形態では、対称性指標マップＭｓ１〜Ｍｓｎの中から、対称性指標Ｓ＝０．９以上となるピクセルが抽出されている。しかし、ピクセルを抽出するための対称性指標Ｓの値は、必ずしも０．９に限られることはなく、必要に応じて、別の値に設定してもよい。

尚、第１の実施形態では、標準偏差σに基づいて、ピーク到達時刻ＴＴＰの値の範囲Ｗ、および、Ｒ_２ ^＊（ｔ）のピーク値Ｖの値の範囲Ｗを決定しているが、標準偏差σ以外の別の値を用いて決定してもよい。

（２）第２の実施形態
第２の実施形態のＭＲＩ装置のハードウェア構成は、図１に示される第１の実施形態のＭＲＩ装置１と同じであるので、ハードウェア構成についての説明は省略する。以下に、第２の実施形態のＭＲＩ装置の動作について、図２２を参照しながら説明する。

図２２は、第２の実施形態のＭＲＩ装置の処理フローを示す図である。

ステップＳ１およびＳ２は、第１の実施形態と同じであるため、説明は省略する。ステップＳ２において信号強度データ系列ＳＤ１（図７参照）を作成した後、ステップＳ３に進む。

ステップＳ３では、信号強度データ系列ＳＤ１に基づいて、造影剤の濃度Ｃｍ（ｔ）の時間変化を表すＣｍデータ系列ＳＤ４を作成する。尚、造影剤の濃度Ｃｍ（ｔ）は、以下の式で表される。
Ｃｍ（ｔ）＝ｋ・ΔＲ_２ ^＊（ｔ） ・・・（１３）
ｋ：定数
ここで、ΔＲ_２ ^＊（ｔ）：式（４）で表されるスピンの緩和速度の差

第１の実施形態で説明したように、ΔＲ_２ ^＊（ｔ）は、信号強度データ系列ＳＤ１が得られていれば、式（４）に従って算出することができる値である。また、定数ｋは、予め決まっている値であり、例えば、ｋ＝２〜３程度の値である。したがって、信号強度データ系列ＳＤ１に基づいて、造影剤の濃度Ｃｍ（ｔ）の時間変化を表すＣｍデータ系列ＳＤ４を作成することができる。信号強度データ系列ＳＤ１は、スライスＡ１〜Ａｎの各領域ごとに求められているので、Ｃｍデータ系列ＳＤ４もスライスＡ１〜Ａｎの各領域ごとに作成される（図２３参照）。

図２３は、作成されたＣｍデータ系列ＳＤ４を示す概念図である。

スライスＡ１〜Ａｎの各領域には、作成されたＣｍデータ系列ＳＤ４が対応付けられている。図２３には、位置座標（ｘ’α、ｙ’α、ｚ’α）における領域Ｒαに対応付けられたＣｍデータ系列ＳＤ４が具体的に示されている。したがって、位置座標が決まれば、対応するＣｍデータ系列ＳＤ４を取り出すことができる。尚、ステップＳ２において既に信号強度データ系列ＳＤ１が対応付けられているので、位置座標ごとに、信号強度データ系列ＳＤ１とＣｍデータ系列ＳＤ４とが対応付けられている。

Ｃｍデータ系列ＳＤ４を作成した後、ステップＳ４に進む。

ステップＳ４では、Ｃｍデータ系列ＳＤ４の特徴量を算出し、特徴量のマップを作成する。第２の実施形態では、Ｃｍデータ系列ＳＤ４の特徴量として、以下の２個の特徴量を算出する。
（１）Ｃｍデータ系列ＳＤ４のピーク部ＰＥにおけるピーク値Ｃmax
（２）ピーク値Ｃmaxにおける時刻（ピーク到達時刻）ＴＴＰ

図２３には、スライスＡｋの領域ＲαにおけるＣｍデータ系列ＳＤ４のピーク値Ｃmax（＝Ｃmaxα）およびピーク到達時刻ＴＴＰ（＝ＴＴＰα）が示されている。しかし、スライスＡｋの他の領域におけるＣｍデータ系列ＳＤ４のピーク値Ｃmaxおよびピーク到達時刻ＴＴＰも算出する。更に、スライスＡｋ以外の他のスライスの各領域におけるＣｍデータ系列ＳＤ４のピーク値Ｃmaxおよびピーク到達時刻ＴＴＰも算出する。

尚、図２３では、Ｃｍデータ系列ＳＤ４の中で、Ｃｍ（ｔ）の最大値をピーク値Ｃmaxαとしている。しかし、Ｃｍデータ系列ＳＤ４に最もよく当てはまるフィッティング曲線を求め、フィッティング曲線の最大値をピーク値Ｃmaxαとするなど、別の方法でピーク値Ｃmaxαを算出してもよい。

図２４は、Ｃｍデータ系列ＳＤ４のピーク値Ｃmaxのマップ、およびピーク到達時刻ＴＴＰのマップを概略的に示す図である。

ピーク到達時刻ＴＴＰのマップ（以下、「ＴＴＰマップ」と呼ぶ）Ｍｔ１〜Ｍｔｎの各ピクセルのピクセル値は、スライスＡ１〜Ａｎの各領域におけるピーク到達時刻ＴＴＰの値を表している。例えば、ＴＴＰマップＭｔｋのピクセルＰαのピクセル値は、スライスＡｋの領域Ｒαにおけるピーク到達時刻ＴＴＰα（図２３参照）を表している。

また、ピーク値Ｃmaxのマップ（以下、「Ｃmaxマップ」と呼ぶ）Ｍｃ１〜Ｍｃｎの各ピクセルのピクセル値は、スライスＡ１〜Ａｎの各領域におけるｔ）のピーク値Ｃmaxの値を表している。例えば、ＣmaxマップＭｃｋのピクセルＰαのピクセル値は、スライスＡｋの領域Ｒαにおけるピーク値Ｃmaxα（図２３参照）を表している。

ＴＴＰマップＭｔ１〜ＭｔｎおよびＣmaxマップＭｃ１〜Ｍｃｎを作成した後、ステップＳ５に進む。

ステップＳ５では、ＴＴＰマップＭｔ１〜Ｍｔｎのピクセルの中から、動脈のピーク到達時刻ＴＴＰを表していると考えられるピクセルを抽出する。ピクセルの抽出方法は、第１の実施形態と同じ方法であるので、説明は省略する。

また、ＣmaxマップＭｃ１〜Ｍｃｎ（図２４参照）のピクセルの中から、動脈における造影剤の濃度Ｃｍ（ｔ）のピーク値Ｃmaxを表していると考えられるピクセルを抽出する。式（１３）に示すように、造影剤の濃度Ｃｍ（ｔ）は、ΔＲ_２ ^＊（ｔ）と比例関係にある。したがって、第１の実施形態において図１３を参照しながら説明した方法と同一の方法で、ピクセルを抽出することができる。

図２４では、ＴＴＰマップＭｔ１〜ＭｔｎおよびＣmaxマップＭｃ１〜Ｍｃｎのピクセルの中で、抽出されたピクセルを、白抜きの四角形で示してある。

ＴＴＰマップＭｔ１〜ＭｔｎおよびＣmaxマップＭｃ１〜Ｍｃｎからピクセルを抽出した後、ステップＳ６に進む。

ステップＳ６では、Ｃｍデータ系列ＳＤ４に基づいて、造影剤の濃度Ｃｍ（ｔ）の２乗値｛Ｃｍ（ｔ）｝^２の時間変化を表す（Ｃｍ）^２データ系列ＳＤ５を作成する（図２５参照）。

図２５は、スライスＡｋの領域Ｒαにおける（Ｃｍ）^２データ系列ＳＤ５を示すグラフである。

図２５には、スライスＡｋの領域ＲαにおけるＣｍデータ系列ＳＤ４（図２３参照）と、Ｃｍ（ｔ）の２乗値｛Ｃｍ（ｔ）｝^２の時間変化を表す（Ｃｍ）^２データ系列ＳＤ５が示されている。

（Ｃｍ）^２データ系列ＳＤ５を作成した後、ステップＳ７に進む。

ステップＳ７では、（Ｃｍ）^２データ系列ＳＤ５のピーク到達時刻ＴＴＰに関する対称性を表す対称性指標Ｓを算出し、対称性指標Ｓのマップを作成する。対称性指標Ｓは、第１の実施形態と同様の方法で時間差ａ１およびａ２を計算することによって算出できる（図２６参照）。

図２６は、対称性指標Ｓのマップを概略的に示す図である。

対称性指標Ｓのマップ（以下、「対称性指標マップ」と呼ぶ）Ｍｓ１〜Ｍｓｎの各ピクセルのピクセル値は、スライスＡ１〜Ａｎの各領域における対称性指標Ｓの値を表している。

対称性指標マップＭｓ１〜Ｍｓｎを作成した後、ステップＳ８に進む。

ステップＳ８では、対称性指標マップＭｓ１〜Ｍｓｎの中から、動脈における対称性指標Ｓを表していると考えられるピクセルを抽出する。動脈における対称性指標Ｓは、一般的に対称性指標Ｓが１に近くなる傾向があるので、第１の実施形態と同様に、対称性指標ＳがＳ＝０．９以上になるピクセルを、動脈の対称性指標Ｓを表していると考えられるピクセルとして抽出する。図２６では、対称性指標マップＭｓ１〜Ｍｓｎの中において、対称性指標ＳがＳ＝０．９以上となるピクセルを、白抜きの四角形で示してある。

ピクセルを抽出した後、ステップＳ９に進む。

ステップＳ９では、先ず、ＴＴＰマップ、Ｃmaxマップ、および対称性指標マップにおいて抽出されたピクセルの中から（図２６参照）、位置座標が共通するピクセルＰα、Ｐβ、およびＰγを検出する。ピクセルＰα、Ｐβ、およびＰγを検出した後、検出されたピクセルの位置座標に対応するＣｍデータ系列ＳＤ４を取り出し、ステップＳ１０に進む。

ステップＳ１０では、取り出されたＣｍデータ系列ＳＤ４を加算平均することにより、動脈入力関数を算出する。

以上説明したように、第２の実施形態では、第１の実施形態とは異なり、造影剤の濃度Ｃｍ（ｔ）を算出し、Ｃｍデータ系列ＳＤ４および（Ｃｍ）^２データ系列ＳＤ５を作成している。しかし、造影剤の濃度Ｃｍ（ｔ）は、式（１３）に示したように、緩和速度の差ΔＲ_２ ^＊（ｔ）と比例関係にある。したがって、ΔＲ_２ ^＊（ｔ）の代わりに、造影剤の濃度Ｃｍ（ｔ）を算出しても、推定誤差の少ない動脈入力関数を得ることができる。

１ ＭＲＩ装置
２ コイルアセンブリ
３ テーブル
５ 受信コイル
６ 制御装置
７ 入力装置
８ 表示装置
９ 被検体
１０ オペレータ
２１ ボア
２２ 超伝導コイル
２３ 勾配コイル
２４ 送信コイル
３１ クレードル
６１ シーケンサ
６２ 送信器
６３ 勾配磁場電源
６４ 受信器
６５ 中央処理装置

Claims (12)

1. 造影剤が注入された被検体の各領域から各時刻における磁気共鳴信号を収集し、収集された前記磁気共鳴信号の信号強度に基づいて前記被検体の血流の動態を解析する血流動態解析装置であって、
   収集された前記磁気共鳴信号に基づいて、前記被検体の各領域ごとに、前記各時刻における磁気共鳴信号の信号強度を算出し、前記信号強度の時間変化を表す第１のデータ系列を作成する第１のデータ系列作成手段と、
   前記第１のデータ系列に基づいて、前記各時刻におけるスピンの緩和速度と、前記造影剤が前記被検体の各領域に到達する前におけるスピンの緩和速度との差を算出し、前記スピンの緩和速度の差の時間変化を表す第２のデータ系列を作成する第２のデータ系列作成手段と、
   前記第２のデータ系列に基づいて、前記各時刻ごとに、前記スピンの緩和速度の差に関連する値を算出し、前記スピンの緩和速度の差に関連する値の時間変化を表す第３のデータ系列を作成する第３のデータ系列作成手段と、
   時間軸方向に関する前記第３のデータ系列の対称性を表す対称性指標を算出する対称性指標算出手段と、
   前記対称性指標に基づいて、前記第２のデータ系列作成手段が作成した第２のデータ系列の中から、動脈入力関数を決定するために使用される前記第２のデータ系列を選択するデータ系列選択手段と、
   を有する血流動態解析装置。
2. 前記第２のデータ系列の特徴量を算出する特徴量算出手段を有する、請求項１に記載の血流動態解析装置。
3. 前記データ系列選択手段は、
   前記特徴量の値と、前記対称性指標の値とに基づいて、前記第２のデータ系列作成手段が作成したデータ系列の中から、動脈入力関数を決定するために使用される前記第２のデータ系列を選択する、請求項２に記載の血流動態解析装置。
4. 前記特徴量算出手段は、
   前記第２のデータ系列のピーク部における第１の時刻を、前記特徴量として算出する、請求項２又は３に記載の血流動態解析装置。
5. 前記対称性指標算出手段は、
   前記第１の時刻に関する前記第３のデータ系列の対称性を表す対称性指標を算出する、請求項４に記載の血流動態解析装置。
6. 前記対称性指標算出手段は、
   前記第３のデータ系列の前記第１の時刻よりも早い時刻に現れるデータ群の重心における第２の時刻と、前記第１の時刻よりも遅い時刻に現れるデータ群の重心における第３の時刻とに基づいて、前記対称性指標を算出する、請求項５に記載の血流動態解析装置。
7. 前記第３のデータ系列作成手段は、
   前記各時刻におけるスピンの緩和速度の差のｐ乗（ｐは２以上の偶数）を、前記緩和速度に関連する値として算出する、請求項１〜６のうちのいずれか一項に記載の血流動態解析装置。
8. 前記第３のデータ系列作成手段は、
   前記各時刻におけるスピンの緩和速度の差のｑ乗（ｑは１以上の奇数）の絶対値を、前記緩和速度に関連する値として算出する、請求項１〜６のうちのいずれか一項に記載の血流動態解析装置。
9. 造影剤が注入された被検体の各領域から各時刻における磁気共鳴信号を収集し、収集された前記磁気共鳴信号の信号強度に基づいて前記被検体の血流の動態を解析する血流動態解析装置であって、
   前記被検体の各領域ごとに、前記各時刻における磁気共鳴信号の信号強度を算出し、前記信号強度の時間変化を表す第１のデータ系列を作成する第１のデータ系列作成手段と、
   前記第１のデータ系列に基づいて、前記各時刻における前記造影剤の濃度を算出し、前記造影剤の濃度の時間変化を表す第２のデータ系列を作成する第２のデータ系列作成手段と、
   前記第２のデータ系列に基づいて、前記各時刻ごとに、前記造影剤の濃度に関連する値を算出し、前記造影剤の濃度に関連する値の時間変化を表す第３のデータ系列を作成する第３のデータ系列作成手段と、
   時間軸方向に関する前記第３のデータ系列の対称性を表す対称性指標を算出する対称性指標算出手段と、
   前記対称性指標に基づいて、前記第２のデータ系列作成手段が作成したデータ系列の中から、動脈入力関数を決定するために使用される前記第２のデータ系列を選択するデータ系列選択手段と、
   を有する血流動態解析装置。
10. 請求項１〜９のうちのいずれか一項に記載の血流動態解析装置を有する磁気共鳴イメージング装置。
11. 造影剤が注入された被検体の各領域から各時刻における磁気共鳴信号を収集し、収集された前記磁気共鳴信号の信号強度に基づいて前記被検体の血流の動態を解析するためのプログラムであって、
    収集された前記磁気共鳴信号に基づいて、前記被検体の各領域ごとに、前記各時刻における磁気共鳴信号の信号強度を算出し、前記信号強度の時間変化を表す第１のデータ系列を作成する第１のデータ系列作成処理と、
    前記第１のデータ系列に基づいて、前記各時刻におけるスピンの緩和速度と、前記造影剤が前記被検体の各領域に到達する前におけるスピンの緩和速度との差を算出し、前記スピンの緩和速度の差の時間変化を表す第２のデータ系列を作成する第２のデータ系列作成処理と、
    前記第２のデータ系列に基づいて、前記各時刻ごとに、前記スピンの緩和速度の差に関連する値を算出し、前記スピンの緩和速度の差に関連する値の時間変化を表す第３のデータ系列を作成する第３のデータ系列作成処理と、
    時間軸方向に関する前記第３のデータ系列の対称性を表す対称性指標を算出する対称性指標算出手段と、
    前記対称性指標に基づいて、前記第２のデータ系列作成処理により作成されたデータ系列の中から、動脈入力関数を決定するために使用される前記第２のデータ系列を選択するデータ系列選択処理と、
    を計算機に実行させるためのプログラム。
12. 造影剤が注入された被検体の各領域から各時刻における磁気共鳴信号を収集し、収集された前記磁気共鳴信号の信号強度に基づいて前記被検体の血流の動態を解析するためのプログラムであって、
    収集された前記磁気共鳴信号に基づいて、前記被検体の各領域ごとに、前記各時刻における磁気共鳴信号の信号強度を算出し、前記信号強度の時間変化を表す第１のデータ系列を作成する第１のデータ系列作成処理と、
    前記第１のデータ系列に基づいて、前記各時刻における前記造影剤の濃度を算出し、前記造影剤の濃度の時間変化を表す第２のデータ系列を作成する第２のデータ系列作成処理と、
    前記第２のデータ系列に基づいて、前記各時刻ごとに、前記造影剤の濃度に関連する値を算出し、前記造影剤の濃度に関連する値の時間変化を表す第３のデータ系列を作成する第３のデータ系列作成処理と、
    時間軸方向に関する前記第３のデータ系列の対称性を表す対称性指標を算出する対称性指標算出処理と、
    前記対称性指標に基づいて、前記第２のデータ系列作成処理により作成されたデータ系列の中から、動脈入力関数を決定するために使用される前記第２のデータ系列を選択するデータ系列選択処理と、
    を計算機に実行させるためのプログラム。

Images