(57)【特許請求の範囲】
【請求項１】
複数の計測量の時系列画像を解析する画像データ解析方法において、
画像抽出手段が、複数の計測量の時系列画像を抽出するステップと、
画像距離設定手段が、前記複数の計測量の時系列画像の各画素において、各画素を固定したときに得られる複数時点での計測量からなるベクトルに距離を設定するステップと、
クラスタリング手段が、前記設定した距離により、前記画素をクラスタリングするステップとからなることを特徴とする画像データ解析方法。
【請求項２】
複数の計測量の時系列画像を解析する画像データ解析方法において、
画像抽出手段が、複数の計測量の時系列画像と、前記複数の画像を標準化するための統一マップを抽出するステップと、
画像変形手段が、前記時系列画像を前記統一マップに基づき変形するステップと、
画像距離設定手段が、前記変形された複数の計測量の時系列画像の各画素において、各画素を固定したときに得られる複数時点での計測量からなるベクトルに距離を設定するステップと、
クラスタリング手段が、前記設定した距離により、前記画素をクラスタリングするステップとからなることを特徴とする画像データ解析方法。
【請求項３】
請求項１に記載の画像データ解析方法において、画像重み付け手段が、前記変形された複数の計測量の時系列画像の各画素において、各画素を固定したときに、時間軸を定義域に複数の計測量を値域とする関数を計算するステップとを有することを特徴とする画像データ解析方法。
【請求項４】
請求項１乃至３いずれかに記載の画像データ解析方法において、画像表示手段が、前記画素のクラスリング結果を表示するステップを有することを特徴とする画像データ解析方法。
【請求項５】
請求項１乃至４いずれかに記載の画像データ解析方法において、フィルタ手段が、前記時系列画像を空間的に平滑化又は時間方向の移動平均の演算のいずれか１つ以上をおこなうことを特徴とする画像データ解析方法。
【請求項６】
時系列画像を解析する画像データ解析システムであって、
複数の計測量の時系列画像と、前記複数の画像を標準化するための統一マップとを抽出する画像抽出手段と、
前記時系列画像を前記統一マップに基づき変形する画像変形手段と、
前記変形された時系列画像の各画素において、各画素を固定したときに得られる複数時点での計測量からなるベクトルに距離を設定する画像距離設定手段と、
前記設定された距離により、前記画素をクラスタリングするクラスタリング手段とを有することを特徴とする画像データ解析システム。
【請求項７】
請求項６に記載の画像データ解析システムにおいて、前記変形された複数の計測量の時系列画像の各画素において、各画素を固定したときに、時間軸を定義域に複数の計測量を値域とする関数を計算する画像重み付け手段を有することを特徴とする画像データ解析システム。
【請求項８】
請求項６又は７に記載の画像データ解析システムにおいて、前記画素のクラスリング結果を表示する画像表示手段を有することを特徴とする画像データ解析システム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、医用画像データの解析方法を搭載した情報システムや医用画像診断装置に関する。
【背景技術】
【０００２】
医用画像診断装置の発展に伴い、様々な物理現象を表す物理量が経時的に得られるようになってきた。特に、Magnetic Resonance Imaging (MRI)では、高精度診断を実現するために、非常に多数の物理量を計測するようになってきている。例えば、物理量としては、核スピン密度、縦緩和時間T1、横緩和時間T2、磁場不均一による見かけ上の横緩和時間T2*、拡散係数、組織・細胞の制限拡散を表す見かけ上の拡散係数（Apparent Diffusion Coefficient、ADC）、血流などの流速、代謝物などの物質の差を表すケミカルシフトなどが挙げられる。これらの組み合わせや、造影剤などの投与による前記物理量の変化など、副次的な情報も含めると、非常に多くの計測量が取得さられるようになってきている。また、従来、主に計測されてきた水素原子核以外の炭素１３やリン３１など他核と呼ばれる原子核について、前記物理量を取得することもできるようになってきた。これらの組み合わせを考えると、得られる計測量は膨大となってきている。また、時間方向については、Echo Planar Imaging (EPI)などデータを高速に取得する方法が開発され、1秒間隔での画像取得も可能となってきている。さらにボリュームデータの高速取得方法の開発とあいまって、時間方向についても、得られる情報は膨大になってきている。
【０００３】
複数の計測量を表す画像を用いた診断技術の例としては、虚血の診断が挙げられる。虚血診断では、核スピン密度、T2、ADC、Perfusion（造影剤を用いた血流量、血液量、平均到達時間）の画像を計測し、組織のダメージを統合的に判断する。他の例としては、がんの確定診断が挙げられる。がんの確定診断では、核スピン密度、T1、T2、ADC、ケミカルシフトによる代謝物濃度などの画像を計測し、悪性度や組織性状などを統合的に判断する。
【０００４】
一方で、単一の計測量ではあるが時系列画像を用いた診断技術の例としては、虚血の診断が挙げられる。虚血診断では、上述した血流量や、血液量、平均到達時間などを計測するために、造影剤投与後のT2画像の変化を追跡している。他の例としては、脳の賦活部位を観測するためのFunctional MRI (fMRI)が挙げられる。ｆMRIでは、感覚刺激や運動刺激とT2*の経時変化との関連性を計測・解析している。
【０００５】
しかし、まだ複数の計測量の時系列画像からの計測・解析技術は確立されていない。一つの理由は、診断に要する労力の飛躍的な増大である。これは、前述したように、観察しなければならない物理量、その変化、一枚の画像のデータ量が増加しており、全てを短時間で見ることは不可能に近くなっているためである。このような診断にかかる労力の増加は、見落としや誤診などを招くこともあるために少しでも軽減することが望まれている。この課題を解決するために、医用画像データを解析する方法がいくつか提案されている。
【０００６】
非特許文献１、２では、脳虚血領域を判別するために、核スピン密度、T2、ADC、Perfusionを用いて画素をクラスタリングする方法が提案されている。また、非特許文献３では、腫瘍の組織判別を行うために、核スピン密度、T2、ADCを用いて画素をクラスタリングする方法が提案されている。画素をクラスタリングする場合、計測された複数の物理量を各軸とする多次元空間に適当な距離を導入し、例えばk-means法などを用いて距離の近い画素を同一な組織としてまとめていく。ただし、これら方法では、時間方向には解析を行わず、一時点での組織判別を繰り返し行っている。
【０００７】
非特許文献４では、アテローム性動脈硬化組織を判別するために、核スピン密度、T1、T2を用いた多次元分析を行う方法が提案されている。本方法では、各画素の物理量だけでなく、各画素の空間的距離も考慮してクラスタリングを行っている。すなわち、非特許文献１、２、３では画素間の空間的な距離は考慮していなかったが、非特許文献４では物理量の軸だけでなく、空間的な距離も軸とした多次元空間に適当な距離を導入している。これにより、離れた領域にある同一性状の組織は別のものとしてクラスタリングされるという効果がある。
【０００８】
非特許文献５では、虚血後のADC変化と乳酸信号量の変化を経時的に計測し、虚血中心やまだ治療が可能と考えられる虚血周辺部を描出できることが示唆されている。これは、虚血領域では、ADCが正常組織の約半分に低下し、乳酸信号が増加していくことを利用している。
【０００９】
【非特許文献１】Carano RAD、 Takano K、 Helmer KG、 Tatlisumak T、 Irie K、 Petruccelli JD、 Fisher M、 Sotak CH. Determination of focal ischemic lesion volume in the rat brain using multispectral analysis. J Magn Reson Imaging 1998;8:1266-1278.
【００１０】
【非特許文献２】Carano RAD、 Li F、 Irie K、 Helmer KG、 Silva MD、 Fisher M、 Sotak CH. Multispectral analysis of the temporal evolution of cerebral ischemia in the rat brain. J Magn Reson Imaging 2000;12:842-858.
【非特許文献３】Carano RAD、 Ross AL、 Ross J、 Williams SP、 Koeppen H、 Schwall RH、 Van Bruggen N、 Quantification of tumor tissue populations by multispectral analysis. Magn Reson Med 2004;51:542-551.
【非特許文献４】Itskovich VV、 Samber DD、 Mani V、 Aguinaldo JGS、 Fallon JT、 Tang CY、 Fuster V、 Fayad ZA. Quantification of human atherosclerotic plaques using spatially enhanced cluster analysis of multicontrast-weighted magnetic resonance images. Magn Reson Med 2004;52:515-523.
【非特許文献５】Takegami T、 Ebisu T、 Bito Y、 Hirata S、 Yamamoto Y、 Tanaka C、 Naruse S、 Mineura K. Mismatch between lactate and the apparent diffusion coefficient of water in progressive focal ischemia. NMR in Biomed 2001;14:5-11.
【発明の開示】
【発明が解決しようとする課題】
【００１１】
前記従来技術では、時系列多値の画像の解析方法が提案されていない。非特許文献５では、虚血後のADC変化と乳酸信号量の変化を経時的に計測し、虚血中心やまだ治療が可能と考えられる虚血周辺部を描出することが示唆はされているが、解析方法は報告されていない。
【００１２】
非特許文献１、２、３、４では、核スピン密度、T1、T2、ADCといった多値画像から組織を判別する方法が提案されているが、時間方向については解析を行っていない。ここで示されている方法では、時間をそれぞれ止めて、各時点での組織判別を行っているのみである。
【００１３】
本発明が解決しようとする課題は、複数の計測量の時系列画像である時系列多値画像を解析し、組織判別を支援する医用画像データ解析方法およびシステムを提供することにある。
【課題を解決するための手段】
【００１４】
本発明は、複数の計測量の時系列画像を解析する画像データ解析方法であって、画像抽出手段が、複数の計測量の時系列画像を抽出するステップと、画像距離設定手段が、複数の計測量の時系列画像の各画素において、各画素を固定したときに得られる複数時点での計測量からなるベクトルに距離を設定するステップと、クラスタリング手段が、設定した距離により、画素をクラスタリングするステップとからなることを特徴とする。
【００１５】
本発明は、複数の計測量の時系列画像を解析する画像データ解析方法であって、画像抽出手段が、複数の計測量の時系列画像と、複数の画像を標準化するための統一マップを抽出するステップと、画像変形手段が、時系列画像を統一マップに基づき変形するステップと、画像距離設定手段が、変形された複数の計測量の時系列画像の各画素において、各画素を固定したときに得られる複数時点での計測量からなるベクトルに距離を設定するステップと、クラスタリング手段が、設定した距離により、画素をクラスタリングするステップとからなることを特徴とする。
また、画像重み付け手段が、変形された複数の計測量の時系列画像の各画素において、各画素を固定したときに、時間軸を定義域に複数の計測量を値域とする関数を計算するステップとを有することを特徴とする。また、画像表示手段が、画素のクラスリング結果を表示するステップを有することを特徴とする。
【００１６】
本発明は、時系列画像を解析する画像データ解析システムであって、複数の計測量の時系列画像と、前記複数の画像を標準化するための統一マップとを抽出する画像抽出手段と、時系列画像を統一マップに基づき変形する画像変形手段と、変形された時系列画像の各画素において、各画素を固定したときに得られる複数時点での計測量からなるベクトルに距離を設定する画像距離設定手段と、設定された距離により、画素をクラスタリングするクラスタリング手段とを有することを特徴とする。
【００１７】
また、変形された複数の計測量の時系列画像の各画素において、各画素を固定したときに、時間軸を定義域に複数の計測量を値域とする関数を計算する画像重み付け手段を有することを特徴とする。また、画素のクラスリング結果を表示する画像表示手段を有することを特徴とする画像データ解析システム。
【発明の効果】
【００１８】
本発明の医用画像データ解析方法によれば、複数の計測値の時間変化から同一な経時変化をとる領域をクラスタとして纏めることが可能となり、全ての画像を詳細に見るための労力を低減することが可能となる。クラスタ結果を観察し、計測値の経時変化をみることで組織判別を支援することが可能となる。また、微小で他とは異なる経時変化をとるような領域を見逃す可能性を低くすることが可能となる。また、ある領域が既知の組織に属する可能性を定量的に評価することも可能となる。すなわち、がんの性状がある物理量の組で表現されている場合、測定された領域がそのがんの性状にどれぐらい近いのかを評価可能となる。
【発明を実施するための最良の形態】
【００１９】
以下、本発明を実施するための最良の形態を説明する。
【実施例１】
【００２０】
本発明の時系列多値画像のデータ解析方法について、図面を用いて説明する。図１に本発明の時系列多値データ解析のフローチャートを表す。
先ず、ステップ１１について説明する。統一マップB(x、y、z)とは、解析する上で空間的な基準を表すものである。x、y、zを空間的な座標を表し、B(x、y、z)は各座標についてある物理量をマップしている関数である。このような統一マップは、複数計測のデータを空間的に統一的に扱うためのものである。例えば、それぞれ異なった空間的な歪がデータに含まれるような計測方法を使用する場合や、異なる患者から標準的な傾向を導き出したい場合に、このような統一マップは必要である。また、このような統一マップの設定は、脳の標準マップのように明示的に示す場合もあれば、一連の計測の内、測定位置を決めるための画像を使用するなど、必ずしも明示的に示さない場合もある。
なお、画素を(x、y、z)で表しているが、２次元空間の(x、y)であっても構わない。
【００２１】
次に、ステップ１２で、設定した統一マップB(x、y、z)にしたがって、時系列多値画像S(m、tm、x、y、z)を変形する。ここで、時系列多値画像とは測定装置を用いてある測定法法により得られる計測データである。x、y、zは空間的な座標を表し、mは計測する物理量、tmは計測する時間を表している。なお、時間tmは絶対的な時刻、もしくは虚血発作開始などの基準となるイベントからの経過時間を表す。計測データ取得に要する計測時間が長い場合には、例えば計測開始時間で代表するなどの操作が為される。
【００２２】
時系列多値画像の変形は、計測方法によって異なるために、各計測量mについてそれぞれ変形を行うことが一般的である。変形には、 （１）空間的な位置の違いを合わせるためのアフィン変換、（２）画素値の個数、解像度の違いを合わせるための補間処理、（３）信号量のむらを合わせるための輝度変換、等を用いる。例えば、（１）については次のような場合がある。B(x、y、z)がSpine Echo (SE)画像であって、S(1、t1、x、y、z)がEcho Planar Imaging (EPI)で取得される画像の場合である。このような場合、EPIでは傾斜磁場の高速なスイッチングによって空間的な位置ずれが生じる。この位置ずれを補正するために、画像重心を計算して平行移動をし、画像内の特徴点を追跡して回転移動と拡大縮小を行うなどの変形処理を行う。例えば、（２）については次のような場合がある。B(x、y、z)がSE画像であって、S(2、t2、x、y、z)がSpectroscopic Imaging (SI)で取得される画像の場合である。このような場合、SIで取得される画像の解像度は低く、画素数も少ない。このため、解像度や画素数を合わせるために線形補間などの空間的な補間処理を行う。例えば、（３）については次のような場合がある。B(x、y、z)が通常のボリュームコイルを用いたSE画像で、S(3、t3、x、y、z)がSNを上げるために局所サーフェスコイルを用いて計測したSI画像の場合である。、サーフェスコイルでは近傍に強い感度を持つので、信号の空間的むらを抑制するために、輝度変換を行う。
【００２３】
なお、抗癌剤による組織変化を解析するような場合、核スピン密度を計測する時系列画像S(1、t1、x、y、z)の変形量を用いて、乳酸を計測する時系列画像S(2、t2、x、y、z)を変形する場合もある。このように、計測量ごとに独立に変形を行えない場合もある。また、この例のように同一の計測量であっても変形量が時間と共に変化する場合もある。
【００２４】
ステップ１３では、変形後の時系列画像S(m、tm、x、y、z)を集めて、各画素(x、y、z)が時系列多値となる時系列多値画像S(m、tm、x、y、z)を作成する。なお、解析プログラムの内部処理としては、S(m、tm、x、y、z)のために新たにメモリを使用することなく、変形したS(m、tm、x、y、z)を用いて処理を行ってもよい。
【００２５】
図２は、この時系列多値画像S(m、tm、x、y、z)の模式図である。時系列多値画像は、関数S(m、tm、x、y、z)で表される。ここで、mは多値画像の各計測量を表し、m=1、2、…Mの値をとる。Mは計測量の個数である。tmは、各計測量mに関する計測時点を表し、計測量mに応じた離散値となる。(x、y、z)は空間的な位置を表す。図面左上のグラフは、横軸にx、縦軸にyをとった空間を表し、各セルは画素を表している。各画素には、計測量であるM個の時系列計測値が割り当てられている。図面では、ある画素(x0、y0、z0)に３種類の時系列計測値が割り当てられている場合を示している。右側一段目はm1の時系列データを示している。S(m1、t、x0、y0、z0)の値が、計測時点t11、t12、t13、t14、t15、t16、t17で構成されている。２段目と３段目は、同一の計測で２種類の計測値を得ている場合に相当する。例えば、Spectroscopic Imaging (SI)画像では様々な代謝物質の信号量を同時に計測することが可能である。このような場合に、乳酸、コリン、クレアチンなどの時系列信号をとっている場合に相当する。
【００２６】
図１のステップ１４は、 (m、tm)の重み付けである。これは、各計測量の重要性を指定するものである。例えば、ADCと乳酸の信号量を二つの計測量とする場合、重み付け無しで比較することは難しい場合もある。例えば、ADCは正常な脳実質であれば約0.8x10-9 m2/sである。乳酸の信号量は任意単位とすると性質も異なり、単位も異なる値を同列で比較することは難しい。これに対して、各計測量に重み付けを行えば、統合的に扱うことが可能となる。もちろん、必ずしもこのステップは必要ではない。例えば、後述するように距離の定義に組み込むことも可能である。
【００２７】
また、各計測量mだけでなく、計測時点tmについても重み付けを設定することができる。例えば、計測時点を細かくとった場合には低く、計測時点を粗くとった場合には高く重み付けを行う。これにより不等間隔で計測を行うような場合に、各計測点の重み付けを適切にし、解析精度を高めることが可能となる。
ステップ１５では、ステップ１４で設定した(m、tm)の重み付けを用いて時系列多値画像S(m、tm、x、y、z)に 重み付け処理を行い、これを時系列多値画像W(m、tm、x、y、z)とする。もちろん、ステップ１４で重み付けが行われない場合には、このステップは必要ない。
【００２８】
ステップ１６では、各画素(x、y、z)を固定した時の時系列多値ベクトルW(m、tm、x、y、z)を、Σtm (m=1、2、…M)次元空間での点と見做し、距離を設定する。図３に本多次元空間の模式図を示す。各矢印は計測量mとその計測時点tmの組を表している。次元数としては全計測時点のΣtm (m=1、2、…M)となる。各画素(x0、y0、z0)を固定したときに相当するベクトルW(m、tm、x0、y0、z0)は、この多次元空間の中の点と考えることができる。この空間に距離を設定する。距離としては、例えばユークリッド距離や、マンハッタン距離、マハラノビス距離を用いることが可能である。
【００２９】
ステップ１７では、設定した多次元空間の距離を用いて、画素(x、y、z)をクラスタリングする。クラスタリング手法としては、k-means、k-nearest neighbor、fuzzy c-clusteringなどの非階層的クラスタリング手法や最近隣法、メジアン法、群平均法などの階層的クラスタリングを利用する。また、組織の性状が大きく異なる場合があるので、このようなクラスタリング手法を数段階に分けて適用してもよい。
【００３０】
最後に、ステップ１８では、画素(x、y、z)のクラスタリング結果を表示する。通常、この表示は統一マップB(x、y、z)を用いて行われる。もちろん、クラスタリング結果を各時系列画像に逆変換して重畳表示するなどしてもよい。また、バッチ処理が行われているような場合には、必ずしも表示は必要なく、解析結果を保存するだけでもよい。
【００３１】
本実施例では、複数の計測量に関する時系列の各計測時点を軸とする多次元空間を考えたが、さらに軸として空間軸x、y、zを入れて距離を設定することも可能である。すなわち、次元数としては3+Σtm (m=1、2、…M)次元の多次元空間で各点をクラスタリングすることとなる。これにより、例えば空間的に近い画素は、多次元空間の中でなるべく近い距離になるように設定することが可能となる。なお、このような空間的に近い画素を同一のクラスタに入れやすくなる方法として、時系列多値画像に空間的に平滑化するフィルタ関数を適用する方法もある。どちらの方法を用いてもよいことは言うまでもない。
【００３２】
図４に、本アルゴリズムを用いた医用画像データ解析方法を搭載した解析システムの機能構成図を示す。本解析システムは大きく三つの手段からなる。一つ目の手段は、本解析プログラムへの画像データの入出力を行う画像データ入出力手段である。これには、画像を蓄積した画像データベースもしくは画像診断装置から画像データを入力する画像抽出手段と、解析した結果の画像を画像データベースに保存する画像保存手段が含まれる。二つ目の手段は、画像をユーザに示す画像表示手段である。これには、統一マップを表示する統一マップ表示手段、解析する元データである時系列多値画像を表示する時系列多値画像表示手段、解析した結果を表示する解析結果表示手段が含まれる。三つ目は画像解析手段である。これには、前記解析アルゴリズムで説明したように、時系列多値画像の統一マップへの変形手段、時系列多値画像の重み付け手段、時系列多値画像の距離設定手段、時系列多値画像の各画素のクラスタリング手段が含まれる。
【００３３】
先ずユーザの入力に従い、画像抽出手段は所定の統一マップと時系列多値画像を画像データベースから抽出する。ここで、統一マップとして、時系列多値画像のうちの一画像、もしくは平均などの計算画像を使う場合がある。このときには、ユーザは直接統一マップの選択を行うことなく、画像抽出手段が所定の方法に従って統一マップを選択もしくは所定の方法で計算する。また、画像診断装置から直接画像を取り込む場合、画像抽出手段が計測画像を抽出する。また、画像データの構造であるが、ＤＩＣＯＭ形式を用いるのが一般的であるが、装置固有のデータ構造であっても構わない。また、時系列多値画像の場合、これを1枚の画像データとするデータ構造で行っても、複数の画像データの集合として扱っても構わない。また、時系列を表すための時間情報は絶対時間であっても、あるイベントからの経過時間であってもよく、そのデータは画像データ構造に含んでも、また別データとして保持していても構わない。ただし、別データとして保持している場合には、それらを読み込む手段が必要となる。
【００３４】
抽出された統一マップと時系列多値画像は、それぞれ統一マップ表示手段と時系列多値画像表示手段によってコンピュータ画面上に表示される。ここで統一マップが明に指定されない場合、特に表示する必要はない。
【００３５】
統一マップと時系列多値画像は画像解析手段に渡される。ここで、ユーザは表示された画像に応じて解析手法やパラメータを選択し、解析を実行する。画像解析手段では必要に応じてユーザの解析手法やパラメータ入力を促すダイアログなどを表示し、前記アルゴリズムにしたがって解析処理を実行する。なお、ここでの解析処理は前記アルゴリズムだけに限ったものではなく、他の実施例で説明するような他の解析処理を行うことも可能である。
【００３６】
結果の画像および時系列グラフなどは解析結果表示手段に渡され、コンピュータ画面上に表示される。ユーザは、必要があれば再度解析を行ったり、解析結果の保存を行ったりする。解析結果の保存が指定された場合、画像保存手段は解析結果を画像データベースに保存する。なお、ここで解析結果の保存先は画像データベースに限らず、印刷やメールや通常のファイルシステムへの保存が行われる場合がある。
【００３７】
図５に、医用画像データ解析方法の画面構成例を示す。本例では画面構成は大きく三つに分かれている。上段は画像解析コントロールパネルで解析対象画像の選択と、解析手法や解析パラメータをユーザが選択できるようにするためのインタフェースである。例えば、左端は統一マップと時系列多値画像を選択するためのボタンが配置されている。その横には解析処理で必要となる重み付け処理の有無や重み付け処理が必要な場合にはその方法が選択できるようにボタンなどが配置される。その右横には時系列多値画像間の距離を設定するためのボタンが配置される。その右横にはクラスタリング手法の選択とそのパラメータを入力するためのボタンと欄が配置される。なお、これらはボタンではなく、メニュー形式や入力欄形式の入力手段であっても構わない。
【００３８】
中段は、ユーザが画像解析コントロールパネルを介して指定したファイルを抽出した画像の表示パネルである。左側に統一マップ、右側に時系列多値画像が表示されている。時系列多値画像では上下で異なる計測量を表し、左右が時間変化を表している。この例では上がＡＤＣで、下が乳酸の信号量を表し、白いほどそれぞれの値が高いことを表している。また、統一マップの輪郭線が重畳表示されている。なお、この例では省力してあるが、これら画像のコントラストなどを制御する設定パネルが配置されてもよい。また、ＡＤＣと乳酸の二値、時間方向が3枚などでなくても良いことは言うまでもない。
【００３９】
下段は解析結果を表示するパネルである。左側は各画素をクラスタリングし、同一クラスタに属する画素は同一色で表示している。右側の折線グラフは、各画素のＡＤＣと乳酸の経時変化を表している。クラスタに従って、折線の色を設定している。
なお、本例は典型的な画面表示を示すもので、このような配置や構成でなくてもよいことは言うまでもない。
【００４０】
図６に、本実施例のアルゴリズムの典型的な変形例を示す。図１と図６との違いはステップ６９にある。それ以外の処理ステップは同等である。
【００４１】
空間(x、y、z)方向や時間tm方向のフィルタリングは、低SN比による解析精度の低下を回避するものである。計測量によってはSN比が低く、このままデータ解析を行うと判別結果がノイズの影響を大きく受けてしまう場合がある。これを抑制するために、空間的に平滑化するフィルタの適用や、時間方向の移動平均の計算、などを行ってもよい。この処理は、図４の画像解析手段にフィルタ手段を設けることよって行う。
【００４２】
時間tm方向の統計値の計算は、解析すべき軸を代表値で置き換え、解析精度を高める場合などに使用する。また、計測量によって計測時点の個数が極端に違う場合などには、計測時点の間引き処理などを行っても良い。また計測時点における計測値の差分などを計算しても良い。これによって経時変化を抽出しやすくなるという効果もある。
【００４３】
さらに、いくつかの計測量mの画像から、計算画像を算出することも可能である。例えば、複数のDiffusion-weighted Imaging (DWI)画像からADC画像を算出する場合のように、指数減衰関数でフィッティングして複数毎の画像を一つに纏めてもよい。この場合、計測量の個数Mや計測時点tmの個数が変化するので、後の処理で次元数を計算する場合に個数を減らすなどの処理が必要である。
【００４４】
なお、ステップ６９で示した処理は、必ずしもフローの中のこの位置で行う必要はない。計測量mで独立に行える空間方向や時間方向の演算処理などは変形前に行うことも可能である。また複数画像の間の画像演算も、それぞれが既に空間的に整合性がおれていれば、変形前に演算処理を行うことも可能である。後で示す実施例４では、一部順序を変更している。
【００４５】
また、ある画素が、複数の計測量の経時変化が既に知られている組織とどれぐらい近いかを定量的に評価することも可能である。この計算には、ステップ１６もしくは６６で設定された距離を使えばよい。
【００４６】
このような解析方法により、通常把握が困難な複数の時系列画像を、単純化できるという利点がある。特に同一な経時変化を見せる画素は同一のクラスタに入るので、代表的な値のみを見ればそれがどのような性状をもつ組織なのか判別するのも容易になるという効果がある。また、完全にこの解析方法だけでは組織判別できない場合でも、この判別結果をもとに元の時系列画像を詳細に観察することも可能になる。さらに、複数の時系列画像を全てみるという多大な労力のもと、見落とす可能性のある微細な他とはことなる経時変化をとる組織などを、抽出しやすくなるという効果もある。
【実施例２】
【００４７】
実施例２は、実施例１のアルゴリズムのうち、ステップ１４以降が変わっている。本アルゴリズムでは、計測時点を多次元空間の軸として各画素をその中の点と考えるのではなく、計測時間から多値への関数空間の中の点と考える。
【００４８】
本アルゴリズムを図７を用いて説明する。
実施例１と異なるステップ７４から説明する。それまでのステップは実施例１、特に図１と同等である。ステップ７４では、先ず各画素(x、y、z)について関数F(m、t)で計測値をフィッティングする。最も単純な場合、関数Fとして、各計測値と計測時点を始点とする階段関数を用いる。ただし、最後の計測の終了時間は、例えば、計測時点の間隔の平均値とするか、計測時点が等間隔であれば、その間隔を用いる。また、関数Fとして、各計測値と計測時点を結ぶ折線グラフで指定される関数としてもよい。また、あるパラメータで規定される関数Fを用いても構わない。例えば、ある計測量に関して、ある時定数で減衰するモデルが適用されると考えられる場合には、計測値を減衰関数でパラメータフィッティングし、その減衰関数をF(m、t)とする。このフィッティングは各計測量mについて独立な場合もあるし、複数の計測量に基づいて行われる場合もある。この選択は適用するモデルに依存する。
【００４９】
次にステップ７５で、図１のステップ１４と同様に計測量の重み付けを設定する。但し、時間tmに関する重み付けは、ここでは必ずしも必要ない。時間方向に重み付けを行いたい場合には、後のステップステップ７７で関数空間の距離を設定する箇所で行うことも可能である。
【００５０】
ステップ７６で、計測量ｍに関する重み付けに従って、多値時間関数画像F(m、t、x、y、z)を重み付けし、多値時間関数画像G(m、t、x、y、z)を計算する。これは、実施例１で述べたように、複数の計測量を統合して扱うためである。この重み付けは、所望のクラスタリング結果が得られるまで、繰り返し変更しながら実施しても良い。
【００５１】
ステップ７７では、各画素(x、y、z)を固定したときの多値時間関数G(m、t、x、y、z)を、時間を定義域にm次元の値を値域とする関数空間の中の点と見做す。この関数空間の中に、距離を定義する。この距離としては、例えば、二乗積分距離や積分距離が利用できる。二乗積分距離とは、二つの関数が与えられたときに、その差分を二乗して時間方向に積分し、積分値の平方根を計算し、それを二つの関数の距離とするものである。
【００５２】
ステップ７８では、設定した関数空間の距離を用いて、画素(x、y、z)をクラスタリングする。クラスタリング手法としては、k-means、k-nearest neighbor、fuzzy c-clusteringなどの非階層的クラスタリング手法や最近隣法、メジアン法、群平均法などの階層的クラスタリングを利用する。また、組織の性状が大きく異なる場合があるので、このようなクラスタリング手法を数段階に分けて適用してもよい。
【００５３】
最後に、ステップ７９で、画素(x、y、z)のクラスタリング結果を表示する。通常、この表示は統一マップB(x、y、z)を用いて行われる。もちろん、クラスタリング結果を各時系列画像に逆変換して重畳表示するなどしてもよい。
【００５４】
このような解析方法により、通常把握が困難な複数の時系列画像を、単純化できるという利点がある。特に同一な経時変化を見せる画素は同一のクラスタに入るので、代表的な値のみを見ればそれがどのような性状をもつ組織なのか判別するのも容易になるという効果がある。また、完全にこの解析方法だけでは組織判別できない場合でも、この判別結果をもとに元の時系列画像を詳細に観察することも可能になる。さらに、複数の時系列画像を全てみるという多大な労力のもと、見落とす可能性のある微細な他とはことなる経時変化をとる組織などを、抽出しやすくなるという効果もある。
【００５５】
実施例１との違いは、関数Fへのフィッティングを使っているため、不等間隔の計測時点でも同様な取り扱いが可能なこと、ノイズに強いことが挙げられる。これにより、計測データによっては、実施例１の方法よりも解析精度を向上できるという効果がある。
【実施例３】
【００５６】
本実施例は実施例１と実施例２で説明した医用画像データ解析方法を搭載するシステムに関するものである。図８に本システムの概略構成図を示す。図８上段は本解析方法が搭載されるデータ解析システム８１を示す。データ解析システムは演算処理を行うCPU、高速に演算処理を行うために一時的にプログラムやデータを記憶するメモリ、プログラムやデータを保存する記憶装置、ユーザからの入力を受け付ける入力装置、ユーザに演算結果などを示す出力装置、医用画像診断装置や他の計算機との通信を行う通信装置から成る。医用画像データ解析方法を記述したプログラムは、記憶装置に保存される。また、医用画像診断装置から通信装置を経由して取得された計測データや計測条件などのデータは、記憶装置に記憶される。ユーザからの本プログラムの起動命令を受けて、本プログラムはメモリに保存される。本プログラムに基づいて、CPUはデータの演算処理を行い、必要であればパラメータ設定などのユーザ入力を促すための入力画面を出力装置に出力し、ユーザからの入力を入力装置で受け付ける。CPUで演算処理された解析結果は出力装置に出力される。
【００５７】
図８下段は、医用画像診断装置８２を表す。ここでは、MRI装置を例に取って説明する。MRIは静磁場発生磁石８３、傾斜磁場コイル８４、RFコイル８５、計算機８６、シンセサイザ、変調器、増幅器、AD変換器で構成される。測定対象８７の核スピンを励起する高周波磁場は、シンセサイザにより発生させた高周波を変調器で波形整形、電力増幅し、RFコイル８５に電流を供給することにより発生させる。傾斜磁場電源から電流を供給された傾斜磁場発生コイル８４は傾斜磁場を発生し、測定対象８７からの磁気共鳴信号を変調する。該変調信号はRFコイル８５により受信され、増幅器で増幅、AD変換器で信号取得された後、計算機８６に入力される。取得されたデータは、計算機でデータ処理され保存される。なお、計算機は予めプログラムされたタイミング、強度で各装置が動作するように制御を行う。保存されたデータは、所定のタイミングで通信装置を介して、データ解析システムに送信される。
【００５８】
なお、本実施例ではデータ解析システム８１と医用画像診断装置８２を別システムとして記述した。もちろん、データ解析システム８１と計算機８６とを同一とすることで、医用画像診断装置自体がデータ解析システムとなることも可能である。
【００５９】
また、本実施例では、MRIを用いて説明したが、X線CTやPET、超音波診断など他の医用診断装置にも適用可能なことは言うまでもない。また、データ解析システムを複数の医用画像診断装置につなぎ、複数の装置で取得されるデータを解析できるようにすることも可能である。この場合、実施例１および２の計測量が複数の装置で得られることになる。
【実施例４】
【００６０】
本実施例では、MRIの計測プロトコルと取得される時系列多値画像、その解析方法および解析結果を示しながら、本発明を説明する。
図９に、本計測プロトコルを示す。上段のグラフは、横軸に時間[分]をとり、計測の流れを示している。計測開始から、まず位置決めなどの計測準備を行い、SE画像を計測する。次に、３０分後からADC計測と、SI計測を交互に３回繰り返す。下にADC計測の撮像シーケンスと、SI計測の撮像シーケンスを示す。シーケンスは横軸に時間、縦軸に各RFやスライス傾斜磁場Gs、リードアウト傾斜磁場Gr、エンコード傾斜磁場Geの強度を表している。ここで、ADC計測は、真ん中のGs、Gr、Geに階段状に示した拡散傾斜磁場を変更しながら複数のDWI計測を行う。ADC画像は、複数のDWI画像から指数関数フィッティングを行って計算される。また、SI計測は、Geの階段状のエンコード傾斜磁場を変更しながら繰り返しデータを計測する。
【００６１】
図８の下段に示したMRI装置は、本シーケンスに記述されたタイミング、強度で稼動し、計測データを取得する。この計測プロトコルによれば、ADC画像とSIで得られる複数の代謝物質画像が、それぞれ３時点の時系列画像として得られる。以降、簡単のため、SI計測で得られる画像を乳酸画像１種類として説明する。すなわち、本計測プロトコルで得られる時系列多値画像は、ADCと乳酸信号の２つの値をとり、時間方向にはそれぞれ３点の計測点をもつこととなる。
【００６２】
なお、本実施例では、本プロトコルにより、脳虚血後のADCと乳酸の経時変化を計測した。
【００６３】
図１０に、本発明を適用した解析フローチャートを示す。比較のため図６のステップも参照しながら説明する。
【００６４】
ステップ１０１は、ステップ６１の統一マップの設定に相当する。統一マップの定義域を、SE画像の最大値の１０％以上となる信号領域に設定する。
ステップ１０２は、複数の計測量の画像から一つの計測量の画像を演算処理するステップ６９に相当する。
【００６５】
ステップ１０３、１０４、１０５は、統一マップへの変形ステップ６８に相当する。SI画像は、画素値が少なく空間分解能も低いため、周波数空間でのゼロフィリングを使って、統一マップに合わせている。
【００６６】
ステップ１０６は、ステップ６３と６６に相当する。今回は、ADCと乳酸の２値で、各計測が３時点からなっているため、６次元の値をとる時系列多値画像となる。
【００６７】
ステップ１０７は、ステップ６４と６５に相当する。ADCは脳実質のほぼ正常値、乳酸はその計測時間での最大値で規格化するように重み付けを行った。ステップ１０８は、ステップ４７のクラスタリングに相当する。ここでは、距離としてユークリッド距離を用い、クラスタリング手法はk-means法を用いた。クラスタリングの個数ｋは、想定される組織の種類である５を用いた。もちろん、距離としては、マハラノビス距離などを用いることも可能であるし、クラスタリングの個数もクラスタを判定しながら、設定することも可能である。
【００６８】
図１１に解析結果を示す。上段は、統一マップで濃淡がクラスタの違いを示している。右側の凡例は、これを見たユーザのデータ解釈を記入したものである。下のグラフは、左にADCの時間推移、右に乳酸の時間推移を表したものである。各折線グラフは、画素１点を表している。凡例にあるユーザの解釈は、このグラフを観察して、記入されたものである。ここでは、正常脳実質と考えられる領域、脳室や脳脊髄液と考えられる領域、虚血中心と考えられる領域、虚血周辺部と考えられる領域、虚血していない領域から血液が流れてくる再灌流が起きていると考えられる領域に分けられる可能性が示されている。
【００６９】
このように、本解析方法によれば、虚血領域を細部まで細かく判別できるという効果がある。特に、再灌流領域や虚血周辺部と考えられる領域の抽出は、複数の時系列画像を眺めているだけでは難しい。本計測プロトコルと本発明の時系列多値画像解析によって、初めて容易になる組織判別方法と考えられる。
【産業上の利用可能性】
【００７０】
医用画像データを解析し、診断の労力を低減することが可能となる。
【図面の簡単な説明】
【００７１】
【図１】本発明の医用画像データ解析方法のアルゴリズム。
【図２】時系列多値画像を表す模式図。
【図３】時系列多値画像の各画素を多次元空間での点と見做した模式図。
【図４】医用画像データ解析方法を搭載した解析システムの機能構成図。
【図５】医用画像データ解析方法を搭載した解析システムの画面表示例。
【図６】本発明の医用画像データ解析方法の別のアルゴリズム。
【図７】本発明の医用画像データ解析方法の別のアルゴリズム。
【図８】本発明の医用画像データ解析方法を搭載するシステムの例。
【図９】時系列多値画像の計測プロトコルの例。
【図１０】本計測プロトコルで計測した時系列多値画像の解析アルゴリズム。
【図１１】本計測プロトコルで計測した時系列多値画像の解析結果。
【符号の説明】
【００７２】
８１・・・データ解析システム、８２・・・医用診断装置、８６・・・計算機。
【図１】
0004622715.tif 000002
【図２】
0004622715.tif 000003
【図３】
0004622715.tif 000004
【図４】
0004622715.tif 000005
【図５】
0004622715.tif 000006
【図６】
0004622715.tif 000007
【図７】
0004622715.tif 000008
【図８】
0004622715.tif 000009
【図９】
0004622715.tif 000010
【図１０】
0004622715.tif 000011
【図１１】
0004622715.tif 000012