JP2008289660A - 脳機能画像分析装置およびその方法並びに脳機能画像分析のためプログラム - Google Patents

Abstract

【課題】時系列の脳機能画像から脳の領域間の関連性を検出する。
【解決手段】本発明の一態様としての脳機能画像分析装置は、時系列の脳機能画像を入力する脳機能画像入力部と、前記脳機能画像における各画素に対応する時系列に基づき前記脳機能画像を複数の領域に分割し、各分割領域を代表する時系列である代表時系列を計算する脳機能画像解析手段と、前記各分割領域の代表時系列をもとに前記各分割領域間の関係性の強さを表す指標を計算し、関係性の強い分割領域間にリンクを設定することにより、各分割領域をそれぞれ１つのノードとしたネットワーク構造を抽出するネットワーク構造抽出手段と、前記ネットワーク構造を出力するネットワーク構造出力部と、を備える。
【選択図】図１

Description

本発明は、脳機能画像分析装置およびその方法並びに脳機能画像分析のためプログラムに関し、たとえばPET（陽電子放出断層撮影法）,fMRI（機能的磁気共鳴法） 等で撮影した脳機能画像の分析に関する。

人間の脳は複数の領域に分かれており、それぞれの領域が何らかの機能を担っていると考えられている。これら複数領域が相互に関連しながら脳全体の機能が作られている。従って脳を理解する上で、個々の領域と相互の関係を理解することは非常に重要である。以下ではこの脳の領域間の関係性をネットワークという言葉で表現する。

従来から脳のネットワークの解明に対して多くの試みがなされてきた。

伝統的には、死んだ脳を解剖し、ニューロンの結合を調べることが多く行われてきた。この方法には物理的な結合の様子を直接見ることが出来るという大きなメリットがある。その反面、脳が活動した状態での動的な関係性をとらえることが出来ないという問題点がある。また、解剖という作業自体に専門的な技能が必要であり、設備と手間を要する。

近年、MRI （磁気共鳴映像法）やX線CT（コンピュータ断層撮影法） といった技術の進歩により、脳構造の詳細な画像がとれるようになってきた。これらの画像から脳の領域を切り出す試みは多く行われている。しかしながら、これらはあくまで静的脳画像に基づくものであり、脳の機能に基づく分割になっていない。

より近年にはfMRI（機能的磁気共鳴映像法）,PET（陽電子放出断層撮影法）,MEG（脳磁図法）,NIRS（近赤外分光法） といった手法により、脳機能を反映した画像が撮影できるようになった。

これらの手法では脳の各部位の活動の様子を表した一連の画像が時系列的に撮影される。ここから脳の部位間の関係性を調べる試みも行われている。文献[FFFDM1997] （R.S.J. Frackowiak, K.J. Friston, C.D. Frith, R.J. Dolan and J.C. Mazziotta, HUMAN BRAIN FUNCTION, Academic Press, 1997. ）にはPATH ANALYSIS による脳機能画像の分割手法が記載されている。この手法は脳の領域間の関係性を調べることを目的としている。しかしながら、脳の複数領域への分割は自動的なものではなく、事前の解剖学的な知識に基づいている。非特許文献１（[DDV2002]）には脳の領域間の関係を因果関係のネットワークとして抽出する方法が示されている。しかしここでも脳を複数領域へ自動的に分割する方法は示されていない。
F.F. Deleus, P. De Maziere and M.M. Van Hulle, Functional Connectivity Modeling in fMRI based on Causal Networks, 2002 IEEE Workshop on Neural Networks for Signal Processing, Martigny, Valais, Switzerland, September 4-6, 2002,pp.505-514.

本発明の目的は、上記の従来技術に鑑み、時系列の脳機能画像から脳の領域間の関連性を自動的に発見する脳機能画像分析装置およびその方法並びに脳機能画像分析のためプログラムを提供する。

本発明の一態様としての脳機能画像分析装置は、
時系列の脳機能画像を入力する脳機能画像入力部と、
前記脳機能画像における各画素に対応する時系列に基づき前記脳機能画像を複数の領域に分割し、各分割領域を代表する時系列である代表時系列を計算する脳機能画像解析手段と、
前記各分割領域の代表時系列をもとに前記各分割領域間の関係性の強さを表す指標を計算し、関係性の強い分割領域間にリンクを設定することにより、各分割領域をそれぞれ１つのノードとしたネットワーク構造を抽出するネットワーク構造抽出手段と、
前記ネットワーク構造を出力するネットワーク構造出力部と、
を備える。

本発明の一態様としての脳機能画像分析方法は、
時系列の脳機能画像を入力する脳機能画像入力ステップと、
前記脳機能画像における各画素に対応する時系列に基づき前記脳機能画像を複数の領域に分割し、各分割領域を代表する時系列である代表時系列を計算する脳機能画像解析ステップと、
前記各分割領域の代表時系列をもとに前記各分割領域間の関係性の強さを表す指標を計算し、関係性の強い分割領域間にリンクを設定することにより、各分割領域をそれぞれ１つのノードとしたネットワーク構造を抽出するネットワーク構造抽出ステップと、
前記ネットワーク構造を出力するネットワーク構造出力ステップと、
を備える。

本発明の一態様としてのプログラムは、
時系列の脳機能画像を入力する脳機能画像入力ステップと、
前記脳機能画像における各画素に対応する時系列に基づき前記脳機能画像を複数の領域に分割し、各分割領域を代表する時系列である代表時系列を計算する脳機能画像解析ステップと、
前記各分割領域の代表時系列をもとに前記各分割領域間の関係性の強さを表す指標を計算し、関係性の強い分割領域間にリンクを設定することにより、各分割領域をそれぞれ１つのノードとしたネットワーク構造を抽出するネットワーク構造抽出ステップと、
前記ネットワーク構造を出力するネットワーク構造出力ステップと、
をコンピュータに実行させる。

本発明により、時系列の脳機能画像から脳の領域間の関連性を自動的に見つけることができる。

（第１の実施形態）
本実施形態は、脳機能画像を機能毎に複数の領域に分割し、それらの相互の関係性を表現したネットワークを抽出することを大きな特徴の１つとする。このような特徴を有する本実施形態としての脳機能画像分析装置の構成を図２に、また本脳機能画像分析装置により行われる処理概要を図１に示す。本脳機能画像分析装置は、脳機能画像を入力する脳機能画像入力部１０と、脳機能画像を機能毎の複数の領域に分割する脳機能画像構造解析手段１１と、それら領域間の関係性を表したネットワーク構造を抽出するネットワーク構造抽出手段１２と、抽出したネットワーク構造を出力するネットワーク構造出力手段１３とを備える。以下脳機能構造解析手段１１およびネットワーク構造抽出手段１２についてそれぞれ詳細に説明する。以下の説明においてたびたび参照される文献の詳細については本欄の最後にまとめて掲載しておく。

まず脳機能画像構造解析手段１１について説明する。

画像分割の手法は数多く提案されている([HG1984])。代表的な手法として、“Region Growing”([Zuck1976]), 境界探査([Davi1975]) などがあげられる。画像分割の対象として脳機能画像を考えた場合、その特徴は次の二つである。
1. 三次元画像であること
2. 一つ一つの画素がそれぞれ一つの時系列になっていること

三次元画像の分割に関しては[LJ1990] に記述がある。これらの文献に示された手法の多くが、上記2の特徴を鑑みた拡張により脳機能画像に応用できる。ここでは一つの実施形態として画素に対応する時系列の類似性に基づく方法を説明する。

図３は、脳機能画像構造解析手段１１の詳細ブロック図である。

画素類似度計算装置２１は脳機能画像の各画素(これに対して一つの時系列が対応する)の間の類似度を計算する。画像分割手段２２は隣接する画素で類似度が高いものを結びつけて、画像を領域分割する。代表時系列抽出手段２３は分割により得られた各分割領域から、その領域に含まれる各画素に対応する時系列を代表する時系列を抽出する。代表時系列モデル２４は代表時系列が満たすような統計的性質を記述したモデルで、代表時系列抽出手段２３による代表時系列抽出の際はこのモデルに則したものが選ばれる。以下、各構成要素をより詳細に説明する。

画素類似度計算装置２１の動作は以下のようになる。既に説明したように、脳機能画像の各画素には一つの時系列(A とする)が対応する。時系列Aの長さをT、時刻t(1≦t≦T)の時系列の値をx_A(t)と表すと、AはT 次元空間上の一つの点

(x_A(1),x_A(2), ・・・,x_A(T))

で表すことが出来る。すべての時系列をこのような点で表現する。このとき、この空間上で近くにある点は時系列として類似していると考えられるから、二つの時系列A、Bの類似度として、この空間上でのユークリッド距離を用いることが出来る。つまり、類似度は以下の式１のように表すことが出来る。また時系列間の類似度の概念を図４に示す。

画像分割手段２２は図５のフローチャートに示すように動作する。まず最終的に作られる領域の数K はあらかじめ指定されているものとする。画像分割手段２２は、画素類似度計算装置２１を用いて画素間の類似度を計算し（Ｓ１１）、計算した画素間の類似度の中で最小のものをとりだし、その二つの画素を統合して１つの領域とする（Ｓ１２）。この操作を領域の数がK になるまで繰り返す（Ｓ１３）。

図６は、画像分割手段２２の動作の具体例を説明するための図である。ここでは最終的な領域の数K は4 とする例を示す。
step1. すべての隣接する画素間の類似度を計算する。初期状態の領域は次のようになる。
{(A), (B), (C), (D), (E), (F), (G), (H), (I)}
step2. 最も類似度が大きいD とG を結合して一つの領域とする。
{(A), (B), (C), (D,G), (E), (F), (H), (I)}
step3. つぎに類似度が大きいD とE を結合して一つの領域とする。
{(A), (B), (C), (D,E,G), (F), (H), (I)}
step4. つぎに類似度が大きいA とD を結合して一つの領域とする。
{(A,D,E,G), (B), (C), (F), (H), (I)}
step5. つぎに類似度が大きいF とI を結合して一つの領域とする。
{(A,D,E,G), (B), (C), (F, I), (H), }
step6. つぎに類似度が大きいA とB を結合して一つの領域とする。
{(A,B,D,E,G), (C), (F, I), (H), }
この段階で四つの領域になったので処理は終了する。

次に代表時系列抽出手段２３の動作を説明する。代表時系列抽出手段２３は、画像分割手段２２により分割された各領域において、代表時系列モデル２４を参照しながら代表時系列を計算する。１つの領域における各画素から１つの代表時系列を計算する様子を図７に示す。

ここでΩがある領域を指すものとして
y_Ω(t) : 領域Ωの代表時系列の時刻t の値
x_i(t) : 領域Ω内のi 番目の画素の時刻t の値
N_Ω: 領域Ω内の画素の数
とする。また、代表時系列モデル２４として各画素の値の代表時系列の値からのずれは偏りの無い正規分布(偏差値σ) に従うと仮定する。つまり

このとき最尤法に従えばy(t)は

できまり、

で与えられる。つまり、最も簡単には代表時系列はその領域内の時系列の平均として求められる。

また別の例では時系列の統計的モデルを参照することも出来る。時系列モデルとしてよく用いられるものとして、MA モデル、ARMA モデル、ARIMA モデルなどが知られている。以下、MA モデルを採用した場合の例を説明する。このモデルは時系列の連続した時刻の間での値の差分が偏りの無い正規分布(偏差値δ) に従うとするものだから

とすればよい。このように複雑なモデルになると、解を与える閉じた数式を書き下すことは出来ないが、数値的な計算は可能である。

次に図２のネットワーク構造抽出手段１２について説明する。ネットワーク構造抽出手段１２の詳細ブロック図を図８として示す。

ネットワーク構造抽出手段１２は代表時系列間の相互の関連性の大きさを表す指標を計算する時系列関連性計算手段３１と、その結果に基づいて二つの領域間にリンクを張るべきかどうかを判断するリンク有無判別手段３２とを備える。

時系列間の関連性としては様々な尺度が利用可能だが、代表的なものをあげると

・時系列間のユークリッド距離
・時系列間の相関(の絶対値)
・時系列間の相互情報量

などがあげられる。
時系列関連性計算手段３１はこのような尺度に基づき代表時系列間の相互の関連性の大きさを表す指標を計算する。

リンク有無判別手段３２は、あらかじめ与えられたしきい値に基づいて、リンクの有無を判断する。つまり、ある二つの代表時系列間の関連性の尺度（指標）がこのしきい値よりも大きければ、そこにリンクを張る。本実施形態において、リンクは方向を有さないものとする。しきい値の決め方として関連性の尺度を直接指定しても良い。

また別の方式では妥当と思われるリンクの数があらかじめ与えるものも考えられる。この場合、関連性の尺度が大きいものから順にリンクを張っていき、与えられたリンク数になったところで停止する。

以上のように、本実施形態によれば、脳機能画像から、脳を機能的なまとまりがある複数の領域に分割することが可能であり、また、その領域の間にある関係性を抽出することが出来る。これによりまず脳機能画像の持つ膨大な情報を分かりやすく整理して示すことが出来る。それと同時に抽出したネットワークの構造を調べることにより、脳の全体的な活動を特徴づけることが出来る。

（第２の実施形態）
第１の実施形態では、まず画像分割を決めた後、分割された各領域の代表時系列を決めるという順に処理が進んだ。しかしながらこの方式では最初の画像分割が代表時系列モデルで指定された時系列の形体と無関係に決定されてしまうので、望ましい代表時系列が得られない可能性がある。そこで別の実施形態として、代表時系列モデルとの兼ね合いを考慮して画像分割を決める方式を説明する。

図９は、本実施形態に係る脳機能画像構造解析手段４１の構成を示すブロック図である。この脳機能画像構造解析手段４１は、図８または後述する図１０のネットワーク構造抽出手段と組み合わせて用いることができる。

本実施形態では、画像の分割が確定してから代表時系列を決める、というような一方方向に事は進まない。代わりに、画像分割手段５２と代表時系列抽出手段５３が協調して画像分割モデル５１と代表時系列モデル５４を同時に最適化するように働く。

ここで画像分割モデル５１とは画像の複数領域への分割の統計的性質を記述したモデルである。今画像がK 個の領域に分割されるものとする。各画素に一つのK 次元ベクトル変数z_iを割り当てて、画素iが領域k (1≦k≦K)に属することを次のように表すものとする。

隣接する画素は同じ領域に属する可能性が高い。このことを確率モデルで表現すると次の式７のようになる。式７は簡単には画像を複数に分割したときにその分割が実現する確率を表している。

ただし積は隣接する画素すべてに対してとるものとする。またZ は分配関数と呼ばれるもので、次の式８で与えら
れる。

和はz₁,・・・,z_N のとり得る値すべてについてとるものとする。Z の計算は実現困難であるが、後の解析においてこの値は必要ではないので、実際にはZ の計算は行わない。

次に領域k の代表時系列が時刻t にとる値をy_k(t)とする。ここで、画素iが時刻t に値x_i(t)をとる確率は次の式９で与えられるものとする。

式８と式９が画像分割モデル５１を構成する。

続いて代表時系列モデル５４について説明する。第１の実施形態のときと同様にここでも時系列モデルとしては一般によく用いられるMA,ARMA,ARIMA が利用可能である。ここでは例としてMA モデルを用いるとすると、次の式１０が代表時系列モデル５４を構成する。

式８、９、１０をあわせることにより、観測値x が与えられたときのy,z の確率モデルは次の式１１で与えられる。

最尤推定法により、y,z の値は

(y, z)=argmax f(y, z|x)・・・式１２

で与えられる。式１２を解いてy（代表時系列の各ｔの値）,z（各画素の属するエリア） の最適値を求めるのは大規模な整数計画問題になり計算量が大きくなる。[ZML1994]では、EM アルゴリズムを用いることにより、妥当な計算量でこの問題を解く手法が示されている。

（第３の実施形態）
図１０は、本実施形態に係るネットワーク構造抽出手段４２の構成を示すブロック図である。本ネットワーク構造抽出手段４２は図３または図９の脳機能画像構造解析手段と組み合わせて用いることが可能である。

時系列記述長計算手段６１はある代表時系列(これをA とする)が与えられたときに、別の代表時系列(これをBとする)の記述長を「関連性」の定量的な尺度として算出するものである。ネットワーク構造を抽出するには、ある領域の代表時系列と別の領域の代表時系列の類似性を調べ、関連性が高いものの間を辺で結ぶというのが基本的な考え方であるが、ここでは「関連性」を、ある代表時系列が与えられたときの別の代表時系列の記述長として定量化している。記述量が正確にどのように定義されるかは後ほど説明する。例えばBの形がA の形に似ていれば、AをもとにBの形を説明するのは簡単にできるであろう。また、BがAと逆位相で変動していれば、やはりAをもとにBの形を説明するのは簡単である。従って、素朴に考えて記述量を「関連性」の定量的な尺度に用いることは自然なことと言える。

AとBの間の記述量はAが与えられたときにBを符号化し、その符号の長さを計算することで求めることが出来る。ランダムな変量を効率よく符号化するには、その確率分布に基づいて行うのが良いことは情報理論ではよく知られた事実である。従ってBを符号化するにはある時刻t におけるAの値a(t)が与えられたときのBの値b(t)の確率分布（時系列間の同時確率密度分布）を求める必要がある。時系列間の同時確率分布密度の算出は同時確率密度分布推定手段６２により行う。

図１１は、横軸にa(t)、縦軸にb(t)をプロットしたものである。これからa(t)とb(t)の同時確率密度分布を求めることが出来る。与えられた点の分布から確率密度分布を求める方法は統計学の分野ではよく研究されている問題である。例えばカーネル法などを用いることが出来る。カーネル法については例えば[Hard1990] に記述されている。

今このようにして求まった同時確率密度分布をp(a,b)としよう。このとき、a が与えられたときのb の確率密度分布は

のように表すことができる。

従ってa(t)が与えられた場合のb(t)の記述量は時系列の各要素ごとに

記述量 =-log{p(b(t)|a(t))} ・・・式１４

で与えられる。

素朴に考えるとこの記述量の計算を一つ一つの時系列要素について積み重ねていけばAが与えられた場合のBの記述量が求まるように思われるが、これは正確ではない。最初に同時確率分布p(a,b)を計算するときに、Bの時系列の値がすべて与えられているとすると、Bの時系列が与えられたもとでBの時系列の記述量を計算するという矛盾したアルゴリズムになってしまうからである。Rissanen はこの問題に対して、予測符号化という解決案を与えている。以下その手法について説明する。

最初の時刻(t =1) においては、a(1) が与えられても、それはb(1) の確率密度分布を考える上で何の助けにもならない。従ってb(1) の符号化に際しては、先験的な分布、例えば一様分布などを用いなければならない。次の時刻(t =2) においては、既に一組の値のペア(a(1),b(1)) が与えられている。ここから、精度は悪くともa とb の同時確率密度分布に関する推定が可能であり、それを用いてa(2) が与えられたときのb(2) の符号長を求めることが出来る。t =1の場合に比較して、符号化の性能があがることが期待できる。以下時刻t が進むにつれて、データが蓄積されa とb の同時確率密度分布の推定の精度が向上するので、符号化の性能もより向上する。このようにしてAが与えられたときのBの符号化を行うのが予測符号化である。a(t),b(t)間の時系列記述長を、時系列記述長計算手段６１と時系列同時確率密度分布推定手段６２との協調により計算する様子を図１２に示しておく。

ここでは一つの時系列Aが与えられときに他の時系列Bの記述長を求める方法を示した。この考え方を拡張して、二つ以上の時系列A₁,A₂,・・・,A_n が与えられたときに、時系列Bの記述長を求めるようにすることも出来る。一般的に参考にする時系列の数n が大きくなるほどBの確率密度分布推定の精度が上がるので、Bの記述長が短くなることが期待できる。

一般にグラフ構造が複雑になる、つまり辺が増える、と個々の代表時系列の記述長は短くなることが期待できる。しかしながら無制限に辺を増やすと、関連性の高さの関係で結ばれる代表時系列の組の数が多すぎて、辺の意味がなくなってしまう。何らかの意味で、辺を増やすことと時系列全体の記述長を減らすことの間に図１３に示すようなトレードオフが無ければならない。ここで最小記述長原理(MDL 原理)と呼ばれる考え方が有効な指針を与えてくれる[Riss1989] 。MDL 原理とは、データを記述する方式として様々な選択肢がある場合、データの記述長を最小にするようなものを採用するべきである、という考え方である。この考え方を今の場合に適用すると、全代表時系列のデータは、グラフ構造と、そのグラフに示された接続関係の元での各時系列の記述で与えられる。すなわち

全データ(ネットワーク全体の)記述長 = グラフ構造の記述長+Σ個々の時系列の記述長
・・・式１５

となる。

グラフ構造の記述長とすべての代表時系列(グラフのすべての頂点（ノード）)の記述長を足しあわせることにより、ネットワーク全体の記述長が求まる。ネットワーク全体の記述長はネットワーク記述長計算手段６３が、上述の時系列記述長計算手段６１と、時系列同時確率密度分布推定手段６２と、以下に説明する非巡回グラフ記述長計算手段６４、非巡回グラフ符号化手段６５とを用いて行う。ネットワーク構造最適化手段６６はネットワーク全体の記述長が最小化するようにグラフ構造を決める。

式１５右辺の二つの記述長のうち、まずグラフ構造の記述長の計算方法を説明する。非巡回グラフ符号化手段６５と非巡回グラフ記述長計算手段６４は協調してネットワークのグラフ構造の記述長を求める。グラフの辺はA → BであればAを使ってBの記述を与えることを示す（後述のように本実施形態では、ネットワークの各ノードのリンクは方向をもつ）。Aを使ってA自身の記述を与えるのは無意味だから、A → B → ・・・ → A のようになることは無い。従ってここで扱うグラフ構造は必ず非巡回(有向)グラフになる。文献[Stein2003] は非巡回グラフを符号化する手法を示している。以下その手法を説明する。

・グラフG は頂点の集合V と、辺の集合E から構成される。
G ≡ (V,E)
頂点はn 個あるものとする。 従って
V = {1,2,・・・,n}
・頂点u の親集合(u を指す辺の出発点となっている頂点の集合)をpa_G(u)とする。
pa_G(u):= {v∈V|(v,u)∈E}

・T(G)はG 中で子が無い頂点の集合を表す。

[非巡回有効グラフのコード化アルゴリズム(by Steinsky)]

非巡回グラフ符号化手段６５は上記のアルゴリズムを用いて、グラフをコード化する。結果として得られるのは

(0,{1},{1,2})

の用な形式の頂点集合の系列である。

非巡回グラフ記述長計算手段６４は次のようにして、この頂点集合の系列の符号長を求める。
・各頂点はlog₂(n)ビットの符号で符号化される。
・各集合に含まれる頂点の数は例えばEllias のω 符号([ 韓・小林1999]) を用いてコード化される。
・系列に含まれる頂点集合の数もやはり例えばEllias のω 符号を用いてコード化される。

非巡回グラフ記述長計算手段６４によってグラフをビット系列にて表現したものを図１４に示す。このビット列から、ビット長を記述長として得ることができる。

このようにして、非巡回グラフ符号化手段６５によって符号化されたグラフ構造は非巡回グラフ記述長計算手段６４によってビット系列に表現され、その記述長が計算できる。

図１０のネットワーク記述長計算手段６３は、代表時系列間のネットワークが後述のネットワーク構造最適化手段６６により与えられると、まず非巡回グラフ記述長計算手段６４を用いてグラフ構造の記述長を求める。次に、グラフ構造の中で親が無い頂点（ノード）に対応する代表時系列の記述長を計算する（たとえば代表時系列における各時刻の値をビット表現する）。更に時系列記述長計算手段６１を用いて、子になっている頂点の記述長を親頂点をもとに順次計算していく。このようにして、ネットワーク全体の記述長が計算できる。

ネットワーク構造最適化手段６６は、ネットワーク全体の記述長を最小化する、代表時系列間のネットワークを求める。本実施形態では、ネットワークの各ノードのリンクは方向をもち、代表時系列Ａから代表時系列Ｂへの方向付き矢印（方向をもったリンク）は、代表時系列Ａを用いて代表時系列Ｂの記述長を求めることを意味している。以下、ネットワーク全体の記述長を最小化する方法を説明する。今の場合、記述長が数式などで与えられているのではないので、グラフ構造に対してどのような操作を行えば記述長を下げられるのかが明確ではない。このような場合に有効な手段としてGA(Genetic Algorithm：遺伝アルゴリズム) が知られている。

一般にGA では探査のための一つの世代が設定される。一つの世代の各メンバは通常、“遺伝コード” と呼ばれるビット列でコード化されている。各メンバの適応度に応じてメンバの遺伝コードを“突然変異”、“交差” などの操作で変更して次の世代を作る。これがGA の標準的な手順である[Gold1989] 。

本例の場合にこのGA を適用した場合、適応度はネットワーク全体の記述長の小ささで測ることが出来る。また、N 個の頂点を持つグラフは、辺が重複しないこと、ノードから自分自身への辺は存在しないこととすれば、存在する可能性のある辺の数はN(N-1)/2 個である。あるグラフに対して実際に存在する辺に対して1、存在しない辺に0を割り振れば、そのグラフはN(N-1)/2 ビットのビット列で表現できる。本例のような非巡回グラフに対してGA の適用を難しくしている点は、あるグラフを表現したビット列に“突然変異”、“交差” の操作を施したとき、非巡回グラフでなくなってしまう可能性があることである。[Novo2003] は“突然変異”、“交差” の代わりに、“influence”,“join” という操作を定義してGA を非巡回有効グラフ構造の最適化に用いる手法を開示している。類似の考え方は[LPYMK1996],[MLD1999],[NK2002] でも示されており、これらの手法はいずれも本発明に適用可能である。

図１５は、ネットワーク構造最適化手段がネットワーク構造を最適化する際の処理の流れを示すフローチャートである。まず、第一世代のグラフ構造をメンバとしていくつか設定し（Ｓ２１）、各メンバの記述長を計算する（Ｓ２２）。そして、各メンバから計算された記述長のうち最小の記述長をl₀に入れる（Ｓ２３）。次に第一世代の各メンバの記述長に応じて次世代のグラフ構造をいくつか設定し（Ｓ２４）、各メンバの記述長を計算する（Ｓ２５）。各メンバから計算された記述長のうち最小のものをl₁に入れ（Ｓ２６）、l₀とl₁との差の絶対値が閾値εより小さいかどうかを判定する（Ｓ２７）。閾値ε以上であるときは、l₁をl₀に入力し（Ｓ２８）、絶対値がεより小さくなるまで、再度Ｓ２４〜Ｓ２６を繰り返す。絶対値がεより小さくなったら処理を終了し、そのときのl₁を得たメンバを採択する。

（第４の実施形態）
先に時系列同時確率密度分布推定手段の実施例として、時系列Bの時刻tにおける値b(t)の確率密度分布を求めるのに、時系列Aの同時刻の値a(t)を用いる方法を示した。これは一番単純な方法であり、様々な拡張が可能である。一つの例として、b(t) の確率密度分布を求めるのに、時系列A の同時刻の値と一つ前の時刻の値a(t), a(t-1) を用いることが考えられる。3 次元空間上の(a(t), a(t-1), b(t)) の分布の例を図１６に示す。この分布を調べることにより、同時確率密度分布p(a, a’, b) を求めることが出来る。このときb の確率分布は次の式で与えられる。

同様な拡張として、b(t) の確率密度分布を求めるのに、時系列A の同時刻の値a(t) と時系列B の一つ前の時刻の値b(t-1) を用いることも考えられる。

なお、本脳機能画像分析装置は、例えば、汎用のコンピュータ装置を基本ハードウェアとして用いることでも実現することが可能である。すなわち、脳機能画像入力手段、脳機能画像構造解析手段、ネットワーク構造抽出手段、ネットワーク構造出力手段、画像分割手段、画素類似度計算装置、代表時系列抽出手段、時系列関連性計算手段、リンク有無判別手段、ネットワーク構造最適化手段、ネットワーク記述長計算手段、時系列記述長計算手段、時系列同時確率密度分布推定手段、非巡回グラフ記述長計算手段、非巡回グラフ符号化手段は、上記のコンピュータ装置に搭載されたプロセッサにプログラムを実行させることにより実現することができる。このとき、本装置は、上記のプログラムをコンピュータ装置にあらかじめインストールすることで実現してもよいし、ＣＤ−ＲＯＭなどの記憶媒体に記憶して、あるいはネットワークを介して上記のプログラムを配布して、このプログラムをコンピュータ装置に適宜インストールすることで実現してもよい。また、代表時系列モデル、画像分割モデルは、上記のコンピュータ装置に内蔵あるいは外付けされたメモリ、ハードディスクもしくはＣＤ−Ｒ、ＣＤ−ＲＷ、ＤＶＤ−ＲＡＭ、ＤＶＤ−Ｒなどの記憶媒体などを適宜利用して記憶することができる。

以下は、本発明の実施の形態の説明において参照する参考文献のリストである。

[FFFDM1997] R.S.J. Frackowiak, K.J. Friston, C.D. Frith, R.J. Dolan and J.C. Mazziotta, HUMAN BRAIN FUNCTION, Academic Press, 1997.

[HG1984] R.M. Haralick and L.G. Shapiro, Image Segmentation Technique, Computer Vision, Graphics, and Image Processing, 29, 100-132, 1985.

[Zuck1976] S.W. Zucker, Region Growing: Childhood and Adolescence, Computer Graphics and Image Processing, 5, 382-399, 1976.

[Davi1975] L.S. Davis, A Survey of Edge Detection Techniques, Computer Graphics and Image Processing 4, pp.248-270,(1975).

[LJ1990] S.P. Liou and R.C. Jain, An Approach to Three-Dimensional Image Segmentation, CVGIP: Image Understanding, Vol. 53,No.3,May,pp237-252,1991.

[ZML1994] J. Zhang, J.W. Modestino and D. A. Langan, “Maximum-Likelyhood Parameter Estimation for Unsupervised Stochastic Model-Based Image Segmantation”, IEEE Trans. on Image Processing,Vol.3,No.4, pp.404-420, 1994.

[Hard1990] W. H¨ardle, Applied Nonparametric Regression, Cambridge University Press, 1990.

[Stein2003] B. Steinsky, Efficient coding of labeled directed acyclic graphs, Soft Computing 7(2003) 350-356.

[DDV2002] F.F. Deleus, P. De Maziere and M.M. Van Hulle, Functional Connectivity Modeling in fMRI based on Causal Networks, 2002 IEEE Workshop on Neural Networks for Signal Processing, Martigny, Valais, Switzerland, September 4-6, 2002,pp.505-514.

[Riss1989] J. Rissanen, Stochastic Complexity in Statistical Inquiry, World Scientific, 1989.

[韓・小林1999] 韓大舜、小林欣吾、「情報と符号化の数理」、培風館、1999.

[Gold1989] D. Goldberg, Genetic Algorithms in Search of Optimization and Machine Learning, Reading:Addison-Wesley.

[Novo2003] A. Novobilski, The Random Selection and Manipulation of Legally Encoded Bayesian Network in Genetic Algorithms, The 2003 International Conference on Artificial Intelligence(ICAI)2003.

[LPYMK1996] P. Larranga, M.Poza, Y. Yurramendi and C.M.H. Kuijpers, Structure Learning of Bayesian NetworksbyGeneticAlgorithms: A Performance Analysis of Control Parameters, IEEE Journal on Pattern Analysis and Machine Intelligence 18(9):912-926.

[MLD1999] J. Myers, K. Laskey and K. DeJong, Learning Bayesian Networks from Incomplete Data Using Evolutionary Algorithms, Proceedings of the Genetic and Evolutionary Computation Conference 1999.

[NK2002] A. Novobilski and F. Kamangar, Genetic Algorithms Based Approach for Discovering Temporal Trends Using Bayesian Networks, The 6^th World Conference on Systematics, Cybernetics and Informatics.

本発明の基本概念を説明する図。 本発明の第１の実施形態としての脳機能画像分析装置の構成を示すブロック図。 図２の脳機能画像分析装置における脳機能画像構造解析手段の詳細構成を示すブロック図。 時系列間の類似度としてのユークリッド距離を説明する図。 図２の脳機能画像分析装置における画像分割手段の動作を説明するフローチャート。 上記画像分割手段の動作の説明図。 代表時系列の説明図。 図２の脳機能画像分析装置におけるネットワーク構造抽出手段の詳細構成を示すブロック図。 第２の実施形態に係る脳機能画像構造解析手段の構成を示すブロック図。 第３の実施形態に係るネットワーク構造抽出手段の構成を示すブロック図。 第３の実施形態に係る時系列間同時確率密度分布推定の説明図。 図１０のネットワーク構造抽出手段における時系列記述長計算手段と時系列同時確率密度分布推定手段による時系列記述長計算処理の流れを示すフロー図。 グラフ構造記述長と全時系列データ記述長のトレードオフの説明図。 非巡回グラフのビット列表現を示す図。 図１０のネットワーク構造抽出手段におけるネットワーク構造最適化手段の計算処理の流れを示すフローチャート。 第４の実施形態に係る時系列間同時確率密度分布推定の説明図。

符号の説明

１０：脳機能画像入力手段
１１、４１：脳機能画像構造解析手段
１２：ネットワーク構造抽出手段
１３：ネットワーク構造出力手段
２１：画素類似度計算装置
２２、５２：画像分割手段
２３、５３：代表時系列抽出手段
２４、５４：代表時系列モデル
３１：時系列関連性計算手段
３２：リンク有無判別手段
５１：画像分割モデル
６１：時系列記述長計算手段
６２：時系列同時確率密度分布推定手段
６３：ネットワーク記述長計算手段
６４：非巡回グラフ記述長計算手段
６５：非巡回グラフ符号化手段
６６：ネットワーク構造最適化手段

Claims (13)

1. 時系列の脳機能画像を入力する脳機能画像入力部と、
   前記脳機能画像における各画素に対応する時系列に基づき前記脳機能画像を複数の領域に分割し、各分割領域を代表する時系列である代表時系列を計算する脳機能画像解析手段と、
   前記各分割領域の代表時系列をもとに前記各分割領域間の関係性の強さを表す指標を計算し、関係性の強い分割領域間にリンクを設定することにより、各分割領域をそれぞれ１つのノードとしたネットワーク構造を抽出するネットワーク構造抽出手段と、
   前記ネットワーク構造を出力するネットワーク構造出力部と、
   を備えた脳機能画像分析装置。
2. 前記脳機能画像解析手段は、隣接する画素間の類似度を、隣接する各画素の時系列から計算し、隣接する画素間の類似度をもとに前記脳機能画像を複数の領域に分割することを特徴とする請求項１に記載の脳機能画像分析装置。
3. 前記脳機能画像解析手段は、前記隣接する画素間の類似度として画素間のユークリッド距離を計算することを特徴とする請求項２に記載の脳機能画像分析装置。
4. 前記脳機能画像解析手段は、最も大きい類似度をもつ画素同士から類似度の降順に順次結合し、あらかじめ指定された領域数に達した時点で結合を停止することにより前記脳機能画像を複数の領域に分割することを特徴とする請求項２または３に記載の脳機能画像分析装置。
5. 前記脳機能画像解析手段は、前記分割領域に含まれる各画素に対応する時系列の平均を前記代表時系列として計算することを特徴とする請求項１に記載の脳機能画像分析装置。
6. 前記脳機能画像解析手段は、前記分割領域に含まれる各画素に対応する時系列と、代表時系列の統計的性質を記述したあらかじめ与えられた代表時系列モデルとから、前記分割領域の代表時系列を計算することを特徴とする請求項１に記載の脳機能画像分析装置。
7. 前記脳機能画像解析手段は、前記脳機能画像の複数領域への分割の統計的性質を記述した画像分割モデルと、代表時系列の統計的性質を記述した代表時系列モデルと、前記脳機能画像における各画素の時系列とに基づいた関数を最適化することにより、前記複数の領域への分割と、前記各分割領域の代表時系列の計算とを同時に行うことを特徴とする請求項１に記載の脳機能画像分析装置。
8. 前記ネットワーク構造抽出手段は、前記分割領域間の関係性の強さを表す指標として、各前記分割領域の代表時系列間のユークリッド距離、各前記分割領域の代表時系列間の相関、または各前記分割領域の代表時系列間の相互情報量を計算することを特徴とする請求項１ないし７のいずれか一項に記載の脳機能画像分析装置。
9. 前記ネットワーク構造抽出手段は、関連性の強さが閾値より大きい領域間にリンクを設定することを特徴とする請求項１ないし８のいずれか一項に記載の脳機能画像分析装置。
10. 前記ネットワーク構造抽出手段は、あらかじめ定められたリンク数に達するまで、関連性の大きい領域間から順にリンクを設定することを特徴とする請求項１ないし８のいずれか一項に記載の脳機能画像分析装置。
11. 前記ネットワーク構造抽出手段は、
    ２つ以上の代表時系列間の同時確率密度分布を推定する同時確率密度分布推定手段と、
    前記同時確率密度分布と、前記２つ以上の代表時系列のうちある代表時系列と異なる他の代表時系列とをもとに前記ある代表時系列の記述長を計算する時系列記述長計算手段と、
    各前記分割領域に対応するノードをもつ非巡回グラフを符号化する非巡回グラフ符号化手段と、
    符号化された非巡回グラフをもとに前記非巡回グラフの記述長を計算する非巡回グラフ記述長計算手段と、
    前記非巡回グラフにおいて親をもたないノードに対応する代表時系列については該代表時系列からその記述長を取得し、親をもつノードに対応する代表時系列については親ノードに対応する代表時系列を前記他の代表時系列とみなして前記時系列記述長計算手段を用いることでその記述長を取得し、前記非巡回グラフにおける各ノードに対応する代表時系列の記述長と、前記非巡回グラフの記述長とを合計したネットワーク記述長を計算するネットワーク記述長計算手段と、
    前記ネットワーク記述長を最小化するように前記非巡回グラフのグラフ構造を最適化するネットワーク構造最適化手段と、
    を備えた請求項１に記載の脳機能画像分析装置。
12. 時系列の脳機能画像を入力する脳機能画像入力ステップと、
    前記脳機能画像における各画素に対応する時系列に基づき前記脳機能画像を複数の領域に分割し、各分割領域を代表する時系列である代表時系列を計算する脳機能画像解析ステップと、
    前記各分割領域の代表時系列をもとに前記各分割領域間の関係性の強さを表す指標を計算し、関係性の強い分割領域間にリンクを設定することにより、各分割領域をそれぞれ１つのノードとしたネットワーク構造を抽出するネットワーク構造抽出ステップと、
    前記ネットワーク構造を出力するネットワーク構造出力ステップと、
    を備えた脳機能画像分析方法。
13. 時系列の脳機能画像を入力する脳機能画像入力ステップと、
    前記脳機能画像における各画素に対応する時系列に基づき前記脳機能画像を複数の領域に分割し、各分割領域を代表する時系列である代表時系列を計算する脳機能画像解析ステップと、
    前記各分割領域の代表時系列をもとに前記各分割領域間の関係性の強さを表す指標を計算し、関係性の強い分割領域間にリンクを設定することにより、各分割領域をそれぞれ１つのノードとしたネットワーク構造を抽出するネットワーク構造抽出ステップと、
    前記ネットワーク構造を出力するネットワーク構造出力ステップと、
    をコンピュータに実行させるための脳機能画像分析のためのプログラム。

Images