(57)【要約】  脳機能データと拡散テンソルデータとを取得し（Ｓ１０）、該拡散テンソルデータに基づき隣接ボクセル間の結合度を算出し（Ｓ３０）、該脳機能データと該隣接ボクセル間の結合度に基づきデータ評価量を構成し（Ｓ４０）、該データ評価量をノンパラメトリック回帰分析し（Ｓ５０）、分析結果に基づき画像を生成し表示する（Ｓ６０，Ｓ７０）。
2007043462.tif 000002
【特許請求の範囲】
【請求項１】
脳の活性部位を特定可能とする脳機能データと前記脳内のプロトンの拡散度を特定可能とする拡散テンソルデータとをボクセル単位で取得する第１のステップと、
前記拡散テンソルデータに基づき隣接ボクセル間の結合度の評価量を構成する第２のステップと、
前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データの解析を行う第３のステップと、
を有することを特徴とする脳機能解析方法。
【請求項２】
請求項１に記載の脳機能解析方法において、前記第３のステップは、前記脳機能データ及びタスクの一方を説明変数とし、前記脳機能データ及び前記タスクの他方を被説明変数とする回帰分析であって、前記脳機能データのボクセル毎の評価量に前記隣接ボクセル間の結合度の評価量を組み込んだ評価量の最小値及び最大値の一方を算出することを特徴とする脳機能解析方法。
【請求項３】
請求項１に記載の脳機能解析方法において、前記第３のステップは、前記脳機能データを検定する手法であって、前記検定の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析方法。
【請求項４】
請求項１に記載の脳機能解析方法において、前記第３のステップは、前記脳機能データを分類する手法であって、前記分類の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析方法。
【請求項５】
請求項１に記載の脳機能解析方法において、前記第３のステップは、前記脳機能データと所定のモデルとの相関をとる手法であって、前記相関の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析方法。
【請求項６】
請求項１に記載の脳機能解析方法において、前記第３のステップは、前記脳機能データから主要な成分を抽出する手法であって、前記抽出の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析方法。
【請求項７】
請求項１に記載の脳機能解析方法において、前記第３のステップにおいて、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データを平滑化することを特徴とする脳機能解析方法。
【請求項８】
請求項１に記載の脳機能解析方法において、前記第３のステップにおいて、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データをクラスタリングすることを特徴とする脳機能解析方法。
【請求項９】
請求項１乃至８のいずれかに記載の脳機能解析方法において、前記第１のステップにおいて取得した拡散テンソルデータに基づき、同ステップにおいて取得した脳機能データから特定される活性部位ボクセルの脳機能データの値を他のボクセルの脳機能データの値に伝播させる前処理を施す第４のステップをさらに備えたことを特徴とする脳機能解析方法。
【請求項１０】
コンピュータを、
脳の活性部位を特定可能とする脳機能データをボクセル単位で取得する脳機能データ取得手段と、
前記脳内のプロトンの拡散度を特定可能とする拡散テンソルデータをボクセル単位で取得する拡散テンソルデータ取得手段と、
前記拡散テンソルデータに基づき隣接ボクセル間の結合度の評価量を構成するデータ評価量構成手段と、
前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データの解析を行うデータ解析手段と、
して機能させることを特徴とする脳機能解析プログラム。
【請求項１１】
請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記脳機能データ及びタスクの一方を説明変数とし、前記脳機能データ及び前記タスクの他方を被説明変数とする回帰分析する手段であって、前記脳機能データのボクセル毎の評価量に前記隣接ボクセル間の結合度の評価量を組み込んだ評価量の最小値及び最大値の一方を算出することを特徴とする脳機能解析プログラム。
【請求項１２】
請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記脳機能データを検定する手段であって、前記検定の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析プログラム。
【請求項１３】
請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記脳機能データを分類する手段であって、前記分類の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析プログラム。
【請求項１４】
請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記脳機能データと所定のモデルとの相関をとる手段であって、前記相関の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析プログラム。
【請求項１５】
請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記脳機能データから主要な成分を抽出する手段であって、前記抽出の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析プログラム。
【請求項１６】
請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析は、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データを平滑化することを特徴とする脳機能解析プログラム。
【請求項１７】
請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データをクラスタリングすることを特徴とする脳機能解析プログラム。
【請求項１８】
請求項１０乃至１７のいずれかに記載の脳機能解析プログラムにおいて、前記コンピュータを、前記第１のステップにおいて取得した拡散テンソルデータに基づき、同ステップにおいて取得した脳機能データから特定される活性部位ボクセルの脳機能データの値を他のボクセルの脳機能データの値に伝播させる前処理を施すデータ前処理手段として機能させることを特徴とする脳機能解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、fMRI（functional Magnetic Resonance Imaging：機能的磁気共鳴撮影）やPET（Positron Emission Tomography：陽電子断層撮影）等の計測手法により得られる各種のデータを用いて脳機能を解析するための脳機能解析方法および脳機能解析プログラムに関する。
【背景技術】
【０００２】
現在、代表的な脳の非侵襲計測方法としては、ｆMRIやPET、MEG（magnetoencephalography：脳磁図）等があるが、このうちデータの空間的分解能が最も高いとされ、広く使用されているのがfMRIである。【０００３】
fMRIは、脳の活性部位を特定可能とする各種の物理量を測定量としてイメージングするものであり、脳機能の計測に対して有効な手法である（非特許文献１及び２を参照）。ｆMRIは、脳の構造をイメージングする解剖MRIの原理と同じく、生体内組織のプロトン密度や縦緩和時間T1、横緩和時間T2を反映するものであるが、特に、脳の活性部位における血流量の増加を捉える点に特徴を有する。脳の活性部位においては局所的に血流量が増加することが知られており、血液中のヘモグロビンは、酸素が結合した状態（酸素化ヘモグロビン）と離れた状態（脱酸素化ヘモグロビン）とでその磁気的性質が異なる。増加した動脈血では磁場を乱す脱酸素化ヘモグロビンの量が少なくなるため、活性部位のfMRI信号（BOLD信号）が増加すると考えられる。したがって、ｆMRIを用いれば、被験者があるタスクを行っているときのBOLD信号の変化を手がかりにして、そのタスクに関係した脳の部位（活性部位）を特定することができる。【０００４】
fMRIにより計測されたBOLD信号の時系列データを解析する際の代表的な手法としては、一般線形モデルに基づいたSPM（Statistical Parametric Mapping）（非特許文献１を参照）や主成分分析、独立成分分析に基づいたデータ解析（非特許文献３を参照）等がある。これらの手法の特徴は、BOLD信号の時系列データを３次元画素（Voxel：ボクセル）毎に個別に統計処理した結果を画像として出力し、脳の活性部位を特定することにある。【０００５】
しかしながら、上記の解析手法においては、データ解析の際に、神経のネットワーク構造が考慮されていないという課題がある。脳は多数の神経細胞がシナプスを介して複雑なネットワークを構成しており、近年の脳科学の知見によれば、脳はこのような神経ネットワークを介して各部位が互いに連携することで全体として一つの高次脳機能を果たしていると考えられている。例えば、被験者があるタスクを行う際に、脳の複数の部位が活性化するといった現象が観測されるが、この現象を上記の解析手法によって解析すると、各活性部位を特定することは可能となるが、活性部位間の結合を特定することは困難である。【０００６】
これにはそもそも以下の理由がある。つまり、上述したようにｆMRIにより計測されるBOLD信号は血流量に基づいているため、脳内において血流量が相対的に多い灰白質（神経細胞の細胞体）の活性は捉えられるが、血流量が相対的に少ない白質（神経細胞の軸索、又は神経線維）の活性はｆMRIでは捉えにくいからである。【０００７】
一方において、神経ネットワーク構造の基盤となるこのような神経線維群の走行方向を捉えるために、近年、MRIによる新たな観測量として生体組織内のプロトンの拡散の度合いを計測するDTI（Diffusion Tensor Imaging：拡散テンソル撮像）法が注目を集めている（非特許文献４を参照）。通常の解剖MRIを利用する場合、神経線維はT1強調画像で高信号、T2強調画像で低信号となる。神経線維がT1強調画像で高信号化する理由はミエリンの存在にあり、ミエリンは二層構造の脂質や巨大タンパク質から成り、神経線維の走行方向に沿って配列するという形態をとる。したがって、神経線維の走行方向ではプロトンの拡散定数が大きく、それに直交する方向では拡散定数が小さいといった異方性が生じる。DTIは、プロトンの拡散を強調するためのMPG（Motion Probing Gradient:傾斜磁場）：
【数１】
2007043462.tif 000003
【０００８】
を印加することによって拡散の異方性を計測する手法であり、例えば、SE（Spin Echo）パルスの前後に拡散検出用のSTG（Stejskal-Tanner Gradient）パルスを加えたST（Stejskal-Tanner）パルス系列によって得られるBOLD信号の強度S’(l,m,k,i)を計測するものである。ここで、l，m，kはボクセルの位置を表す正の離散変数であり、それぞれ、ボクセルのX方向、Y方向、Z方向の座標を表す。また、iは、測定時間を表す正の離散変数である。該BOLD信号の強度S’(l,m,k,i)は、
【数２】
2007043462.tif 000004
【０００９】
と書くことができ、拡散強調画像を生成する際には式（１）における拡散テンソル
【数３】
2007043462.tif 000005
【００１０】
がデータ解析の対象となる。ここで、ρ’(l,m,k,i)はMPGを印加しない場合のBOLD信号(通常の脳機能解析におけるデータ解析の対象)の強度を表し、bはMPGの強さを表すパラメータである。なお、ρ’(l,m,k,i)は、
【数４】
2007043462.tif 000006
【００１１】
と表される。ここで、f(ν)は流速、TRは繰り返し時間、TEはエコー時間、そしてξ’(l,m,k,i)はプロトン密度を表す。【００１２】
現在、fMRIにより計測されて得られたBOLD信号の時系列データρ’(l,m,k,i)を解析した後に、脳の関心領域（Region of Interest:ROI）を拡散テンソルデータD(l,m,k)によって接続するといった脳機能解析方法が用いられつつある（非特許文献５を参照）。
【非特許文献１】“Human Brain Function：2nd-Ed.”, Richard S. J. Frackowiak, et al, ELSEVIER ACADEMIC PRESS, 2004
【非特許文献２】“Image of Mind”, M. I. Posner and M. E. Raichle, W H Freeman & Co, 1997
【非特許文献３】“Independent Component Analysis: Theory and Applications”, T.-W. Lee, Kluwer Acadmic, 1988
【非特許文献４】“これでわかる拡散MRI”青木 茂樹, 阿部 修，秀潤社，2002
【非特許文献５】“Combined functional MRI and tractography to demonstrate the connectivity of the human primary motor cortex in vivo”， Guye M，et al.， Neuroimage， Vol.19， pp.1349-1360， 2003
【発明の開示】
【００１３】
しかしながら、上記の脳機能解析方法においては、BOLD信号の時系列データρ’(l,m,k,i)自体の解析には拡散テンソルデータD(l,m,k)が用いられていないため、その有効性には疑問が残る。
【００１４】
本発明は、上記の課題を鑑み為されたものであり、その目的とするところは、ｆMRIやPET等の非侵襲計測方法を用いて計測された脳機能データと神経線維群の走行方向を特定可能とする拡散テンソルデータとに基づき脳の活性部位間の結合構造をも考慮した脳機能データの解析を可能とする脳機能解析方法および脳機能解析プログラムを提供することにある。
【００１５】
本発明の一つの側面は、脳の活性部位を特定可能とする脳機能データと前記脳内のプロトンの拡散度を特定可能とする拡散テンソルデータとをボクセル単位で取得する第１のステップと、前記拡散テンソルデータに基づき隣接ボクセル間の結合度の評価量を構成する第２のステップと、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データの解析を行う第３のステップと、を有することを特徴とする脳機能解析方法を提供することである。
【００１６】
本発明の他の側面は、コンピュータを、脳の活性部位を特定可能とする脳機能データをボクセル単位で取得する脳機能データ取得手段と、前記脳内のプロトンの拡散度を特定可能とする拡散テンソルデータをボクセル単位で取得する拡散テンソルデータ取得手段と、前記拡散テンソルデータに基づき隣接ボクセル間の結合度の評価量を構成するデータ評価量構成手段と、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データの解析を行うデータ解析手段と、して機能させることを特徴とする脳機能解析プログラムを提供することである。
【図面の簡単な説明】
【００１７】
【図１】図１は、本発明の第１の実施形態に係る脳機能解析装置の概略構成を示したブロック図である。
【図２】図２は、典型的なｆMRI計測の方法を例示した図である。
【図３】図３は、典型的なｆMRI画像（２次元スライス画像）を例示した図であり、図３Ａは寝台上の被験者の頭部を表し、図３Ｂは頭部の２次元断面（２次元スライス画像）を構成するボクセルを表す。
【図４】図４は、２次元スライス画像を構成するボクセル上の脳機能情報の時系列データと拡散テンソルデータを例示した表である。
【図５】図５は、図４に示した脳機能情報の時系列データの一例を示した表であり、図５Ａは実測値を表し、図５Ｂは実測値を閾値に基づき二値化した値を表す。
【図６】図６は、図４に示した脳機能情報の時系列データと拡散テンソルデータとを立体構成した図であり、図６Ａは寝台上の被験者の頭部を表し、図６Ｂは図４に示した脳機能情報の時系列データと拡散テンソルデータとからｋ番目の２次元スライス画像Skを生成する手順を表す。
【図７】図７は、図６に示した２次元スライス画像Skから３次元構成された頭部画像を例示した図であり、図７Ａは寝台上の被験者の頭部を表し、図７Ｂは２次元スライス画像Skを集めて３次元的に構成された頭部画像を表す。
【図８】図８は、図１に示した脳機能解析装置によって行われるデータ解析の処理手順の一例を示したフローチャートである。
【図９】図９は、本発明に係るデータ前処理の手法の一例を示した図であり、図９Ａはあるボクセルに隣接する２つのボクセルが共に灰白質と見なされる領域におけるデータ前処理を表し、図９Ｂはこれら２つのボクセルのうちの１つだけが灰白質と見なされる領域におけるデータ前処理を表し、図９Ｃはこれら２つのボクセルが共に灰白質とは見なされない領域におけるデータ前処理を表す。
【図１０】図１０は、本発明の第２の実施形態に係る脳機能解析装置の概略構成を示したブロック図である。
【図１１】図１１は、図１０に示した脳機能解析装置によって行われるデータ解析の処理手順の一例を示したフローチャートである。
【図１２】図１２は、図１０に示したデータ平滑化手段によって行われる平滑化の手法を例示した図である。
【図１３】図１３は、第２の実施形態の変形例に係る脳機能解析装置の概略構成を示したブロック図である。
【図１４】図１４は、図１３に示した脳機能解析装置によって行われるデータ解析の処理手順の一例を示したフローチャートである。
【図１５】図１５は、本発明の第３の実施形態に係る脳機能解析装置の概略構成を示したブロック図である。
【図１６】図１６は、図１５に示した脳機能解析装置によって行われるデータ解析の処理手順の一例を示したフローチャートである。
【図１７】図１７は、図１５に示したクラスタリング手段３００によって行われるクラスタリングの手法を例示した図であり、図１７Ａはクラスタリング処理を行う前の結合度ベクトルの値を表し、図１７Ｂはクラスタリング処理を行った後の結合ベクトルの値を表す。
【図１８】図１８は、第３の実施形態の変形例に係る脳機能解析装置の概略構成を示したブロック図である。
【図１９】図１９は、図１８に示した脳機能解析装置によって行われるデータ解析の処理手順の一例を示したフローチャートである。
【図２０】図２０は、本発明の第４の実施形態に係る脳機能解析装置の概略構成を示したブロック図である。
【図２１】図２１は、図２０に示した脳機能解析装置によって行われるデータ解析の処理手順の一例を示したフローチャートである。
【図２２】図２２は、第４の実施形態の変形例に係る脳機能解析装置の概略構成を示したブロック図である。
【図２３】図２３は、図２２に示した脳機能解析装置によって行われるデータ解析の処理手順の一例を示したフローチャ−トである。
【図２４】図２４は、被験者に聴覚刺激を用い単純復唱を行わせた場合の各種脳機能解析を行った結果を示した図であり、図２４Ａ，２４Ｂは共にSPMを用いたT検定（図２４ＡはT検定の閾値を修正した場合，図２４ＢはT検定の閾値を無修正の場合）の結果を表し、図２４ＣはSPMとトラクトグラフィーとを組み合わせた脳機能解析手法を用いた解析結果を表し、図２４Ｄは本発明の第１の実施形態に係る脳機能解析方法を用いた解析結果を表す。
【発明を実施するための最良の形態】
【００１８】
本発明は、脳の活性部位を特定するために、MRI等の非侵襲計測装置から脳機能データと拡散テンソルデータとをボクセル単位で取得し、該拡散テンソルデータから隣接ボクセル間の結合度の評価量を構成し、該隣接ボクセル間の結合度の評価量を用いて、該脳機能データの解析を行うことを特徴とするものである。上記脳機能データの隣接ボクセル間の結合度の評価量はデータ解析の手法に依存して様々考えられるし、また該隣接ボクセル間の結合度の評価量をどのように用いるかについてもデータ解析の手法に応じて様々考えられる。【００１９】
そこで、以下に、本発明の実施形態を、図面を用いて詳細に説明する。【００２０】
［第１の実施形態］
図１は、本発明の第１の実施形態に係る脳機能解析装置の概略構成を示したブロック図である。第１の実施形態の脳機能解析装置１０は、MRI装置５０で計測された脳機能情報の原時系列データρ’(l,m,k,i)（例えば、(3)式）を取得する脳機能データ取得手段１と、同じくMRI装置５０で計測された拡散テンソルデータD(l,m,k)（例えば、(2)式）を取得する拡散テンソルデータ取得手段２と、脳機能データ取得手段１で取得した脳機能情報の原時系列データρ’(l,m,k,i)に対して前処理を施すデータ前処理手段３と、拡散テンソルデータ取得手段２で取得された拡散テンソルデータD(l,m,k)から隣接ボクセル間の結合度を表す結合度ベクトル
【数５】
2007043462.tif 000007
【００２１】
を算出するボクセル間結合度算出手段４と、データ前処理手段３で前処理された脳機能情報の時系列データρ(l,m,k,i)とボクセル間結合度算出手段４で算出された結合度ベクトル
【数６】
2007043462.tif 000008
【００２２】
とから所定の評価量Qを構成するデータ評価量構成手段５と、該データ評価量Qの極値（本実施形態では最小値）を求めるための計算を行うデータ解析手段６と、データ解析手段６で算出されたデータに基づき画像を生成する画像生成手段７と、画像生成手段７によって生成された画像を表示する画像表示手段８と、上記各手段２〜８で取得、算出、構成、解析、及び生成された各種データを記録するサブメモリと、以降で詳述する各ステップを実行するためのコンピュータに読み取り可能な脳機能解析用プログラムを記憶するメインメモリとから成る記憶手段９と、記憶手段９から読み出された該脳機能解析用プログラムに従って上記各手段２〜９を制御するCPU(Central Processing Unit：中央演算処理装置)等を備える。なお、上記の各記号の意味は後に説明する。【００２３】
ここで、MRI装置５０は、核磁気共鳴（Nuclear Magnetic Resonance:NMR）を起こすための静磁場を生成する静磁場コイル、共鳴周波数の高周波を照射し、共鳴信号を検出する高周波コイル、該共鳴信号に位置情報をエンコードするための傾斜磁場を生成する傾斜磁場コイル等から成るマグネットアセンブリと、これらのコイルの通電を制御するシステムコントローラ等から構成されるものであり、操作者の要求に応じて各種fMRIデータ（血流量に関連する脳機能情報の原時系列データρ’(l,m,k,i)や神経繊維の走行方向に関連する拡散テンソルデータD(l,m,k)等）を生成し、これらのデータを脳機能解析装置１０に送信するものである。また、本発明においては、脳機能データの原時系列データρ’(l,m,k,i)を必ずしもMRI装置５０から取得する必要はなく、その代わりに、同じく各種の脳機能データを生成可能なPET等の非侵襲計測装置から取得する構成としてもよい。【００２４】
また、画像表示手段８は、例えば、CRTディスプレイ、TFTディスプレイ、プラズマディスプレイなどの各種ディスプレイや、インクジェットプリンタ、レーザープリンタなどの各種プリンタなどが使用可能である。また、記憶手段９は、例えば、RAM（Random Access Memory）やROM（Read Only Memory）等から構成されている。さらに、記録手段９のサブメモリとメインメモリとを別体として構成し、メインメモリ部分を磁気ハードディスク、フロッピー(登録商標)ディスク、CD-ROMなどの光ディスク、磁気テープ、メモリチップ等に記憶させてもよい。【００２５】
本実施形態においては、脳機能解析装置１０を、本実施形態の主要構成手段１０Ａと画像表示手段８及び記憶手段９とを一体型としたディスプレイコンソール形式としたが、画像表示手段８、或いは、画像表示手段８及び記憶手段９をそれぞれ独立した画像表示装置及び記憶装置として主要構成手段１０Ａから切り離した構成としてもよい。いずれの構成においても、脳機能解析装置１０はコンピュータによって実現されるものであり、上記各手段２〜９は、CPU（Central Processing Unit：中央演算処理装置）１００により記憶手段９から読み出された脳機能解析用プログラムに従って制御される。【００２６】
ここで、コンピュータとは、構造化された入力を所定の規則に従って処理し、処理した結果を構造化して出力する装置のことを指し、例えば、汎用コンピュータ、スーパーコンピュータ、メインフレーム、ワークステーション、マイクロコンピュータ、サーバ等が含まれる。また、通信ネットワーク（例えば、イントラネット、ローカルエリアネットワーク（LAN）、ワイドエリアネットワーク（WAN）、及びこれらの組み合わせから成る通信ネットワーク）を介して接続された２つ以上のコンピュータから成る構成（例えば、分散コンピュータシステム）であってもよい。【００２７】
次に、本発明の理解を容易にするために、図２乃至図７を用いて、本発明による脳機能解析の流れの概略を説明する。図２は、典型的なｆMRI計測の方法を例示したものであり、図３は、典型的なｆMRI画像（２次元スライス画像）を例示したものであり、図３Ａは寝台上の被験者の頭部を表し、図３Ｂは頭部の２次元断面（２次元スライス画像）を構成するボクセルを表す。【００２８】
例えば、被験者にあるタスク（例えば、指のタッピング）を行わせた際に、そのタスクに関連すると考えられる脳の部位をｆMRI計測によって特定することを考える。図２に示す例では、１回の計測において、一定時間Tのタスクと同じ一定時間Tの休憩（レスト）とを１セットとして３回繰り返す場合を示している。ここで、横軸は時間軸であり、“ON”はタスクを表し、“OFF”はレストを表す。通常、ｆMRI計測は１回の実験において数十回行われる。この例では、ｆMRI計測は２４回行われるわけだが、タスクの立上がりと立下りの画像は計測時刻のずれを考慮して通常使用しないので、実際には、有意なｆMRI計測は１８回（図中の１８本の太線に対応）行われたことになる。なお、ti（for i=1,2,…,24）は計測の実時間を離散的に示している。以降では、tiを単に“i（for i=1,2,…,I）”（正の離散整数）と表す。【００２９】
このようなｆMRI計測の結果、図３Ｂに示すような２次元スライス画像Skが得られる。ここで、kは１回の実験によって得られる２次元スライス画像の枚数を表し、通常のｆMRI計測においては２０〜３０枚の２次元スライス画像が得られる（図中、K=20〜30）。これらの２次元スライス画像Sk（for k=1, 2,…,K）は、図３Ｂにおいては平面的に記載されているが、実際にはボクセルと呼ばれる３次元の立体画素から構成されている。同図では、一例として６４×６４個のボクセルから成る２次元スライス画像を示したが、もちろん、その他の個数のボクセルに分割した２次元スライス画像であってもよい。一般に、任意の２次元スライス画像Sk上のボクセルは２つの正の離散変数l,mを用いて(l,m)と表すことができるが、所定の規則で番号付けすることによってこれを単一の正の離散変数jで表すこともできる。つまり、数学的には、jと(l,m)とを該所定の規則によって１対１対応（j≡(l,m)（for Sk））させることが可能である。この例では、lの上限L及びmの上限Mが共に６４（L=64,M=64）であることから、j=1,2,…,4096となる。以降においては、これらの表記を適宜使い分けることにする。また、図３Ｂの丸は頭部の断面形状の輪郭を表している。【００３０】
図４は、脳機能データ取得手段１で取得し、データ前処理部３で前処理（後述）を施した、ある２次元スライス画像Skを構成するボクセルj上の脳機能情報の時系列データρ(j,k,i)（≡ρk(j,i)）と拡散テンソルデータ取得手段２で取得されたD (j,k)（≡Dk(j)）を例示したものである。背景技術で述べたように、DTI計測によって拡散強調画像を生成する際には、典型的にはSEパルスの前後に拡散検出用のSTGパルスを加えたSTパルス系列でイメージングした結果として、該２次元スライス画像Sk上の脳機能情報の原時系列データρ’(l,m,k,i)（≡ρ’k(l,m,i)）と拡散テンソルデータD(l,m,k) （≡Dk(l,m)）とが得られる。ρk(j,i)は、この原時系列データρ’k(l,m,i)に対して後述する所定の前処理を施した結果得られる脳機能情報の時系列データを表す。つまり、図中のρk(j,i)は、ある２次元スライス画像Skのボクセルjにおける時刻iの脳機能情報を表す。また、Dk(j)は、ある２次元スライス画像Skのボクセルjにおける拡散テンソル情報を表す。ρk(j,i)は具体的にはスカラーとして表され、また、Dk(j)は式（２）に示したように行列で表現され、例えば、
【数７】
2007043462.tif 000009
【００３１】
のような値をとる。この例で示すように、拡散テンソルは通常、対称テンソル（非対角成分が対称）なので、実質的には６個の成分を有することになる。また、図中のCLASSとは、タスクを行っているか否かを表したものであり、“ON（=1）”及び“OFF（=0）”はそれぞれ図２の“ON”及び“OFF”に対応する。【００３２】
図５は、図４に示した脳機能情報の時系列データρk(j,i)の一例を示した表である。図５（Ａ）に例示するように、実際の脳機能情報の時系列データρk(j,i)は連続値を取るが、例えば図５（Ｂ）に示すように、値“4000”を活性の閾値として設定して、それよりも大きい値を取るボクセル値は活性“A”、小さい値を取るボクセル値は不活性“I”というように前処理の段階でさらに二値化してもよい。【００３３】
図６は、図４に示した拡散テンソルデータDk(j)を用いて、同じく図４で示した脳機能情報の時系列データρk(j,i)とを立体構成した図であり、図６Ａは寝台上の被験者の頭部を表し、図６Ｂは図４に示した脳機能情報の時系列データと拡散テンソルデータとからｋ番目の２次元スライス画像Skを生成する手順を表す。図６Ｂのボクセルj1及びj2上の５角形は、図４に示した脳機能情報の時系列データρk(j,i)及び拡散テンソルデータDk(j)が付随した各ボクセル上の内部空間（つまり、ボクセル上のスカラー場及びテンソル場）を例示したものである。【００３４】
例えば、背景技術で述べたSPMにおいては、MRI装置５０から取得した脳機能情報の原時系列データρ’k(l,m,i)に対して、
（ａ）整列（Realignment）：先頭の２次元スライス画像（例えば、S1）に、以降の画像をそろえることで、測定中の頭の動きに伴う位置補正をし、動きによる擬似信号を取り除き、
（ｂ）空間的正規化（Spatial Normalization）：複数の被験者のデータをまとめたり比較したりするために、各被験者のデータをTalairach標準脳に合わせ、
（ｃ）スムージング（Smoothing）：ノイズが多く含まれる原時系列データにガウス型フィルターを適用することで、空間的分解能を下げることなく、解析の感度を向上させる、
といった前処理を施した後に、ボクセル毎に各種の検定等を行い、タスクに伴う有意な活性部位を、被験者の個人データで、さらには複数の被験者から成るグループ間で検定する。さらに、SPMでは、活動推定のモデルを事前に用意し、上記のようにして前処理された機能情報の時系列データρk (j,i)（=ρk (l,m,i)）がどれだけモデルに照合するかを、拡散テンソルデータDk(l,m) を用いることなく、一般線形モデルを使いパラメータ推定する。【００３５】
一方、本発明では、脳機能情報の原時系列データρ’k(l,m,i)に対して、拡散テンソルデータ取得手段２によって取得された拡散テンソルデータDk(j)（=Dk(l,m)）を用いて、上記の前処理以外にさらに以下で詳述するような前処理（ｄ）を施す。そして、該拡散テンソルデータDk(j)（=Dk(l,m)）を用いて、前処理された脳機能情報の時系列データρk(j,i)（=ρk(l,m,i)）に以下で詳述するようなデータ処理を施すことによって画像表示の対象となる２次元画像データak(j)（=ak(l,m)）を生成する。【００３６】
なお、上記では簡単のために２次元スライス画像を例にとり説明したが、実際には、図７Ａ，７Ｂに示すような頭部の立体画像データak(j)（=a(j,k)=a(n)）を３次元的にデータ解析する。その際にも、以降で詳述するような、拡散テンソルデータDk(j)を用いたデータ処理が施される。ここで、“n”とは、立体構成された全ボクセルに対して所定の規則で番号付けしたボクセルの位置を表す正の離散変数である。図中では、XY座標上の２次元スライス画像の位置を表す変数jと３次元目の方向（z方向）を表す変数kとを分けて記載したが、これを上記した１変数“n”で表すこととは該所定の規則の下で数学的に同値である（n≡(j,k)）。なお、以降においては、同一の３次元ボクセルを表すのに１変数“n”を用いた表記とXYZ座標系における３変数(j,l,k)を用いた表記とを適宜使い分ける（n≡(l,m,k)≡(j,k)）。【００３７】
以上の前提の下、図１に示した脳機能解析装置１０によって行われるデータ解析方法について詳細に説明する。図８は、図１に示した脳機能解析装置１０によって行われるデータ解析の処理手順を示したフローチャートである。【００３８】
ステップＳ１０において、脳機能データ取得手段１は、MRI装置５０で計測された脳機能情報の原時系列データρ’(l,m,k,i)を取得する。同ステップにおいて、拡散テンソル情報取得手段２は、同じくMRI装置５０で計測された拡散テンソルデータD(l,m,k)を取得する。【００３９】
次に、ステップＳ２０において、データ前処理手段３は、上記のようにして取得した脳機能情報の原時系列データρ’(l,m,k,i)に対して、SPMの上記手法（ａ）〜（ｃ）に代表されるような位置補正、ノイズ除去等の前処理を施す。【００４０】
同ステップにおいて、データ前処理手段３は、上記のようにして通常の前処理が施された脳機能情報の時系列データρ’(l,m,k,i)に対してさらに前処理（ｄ）（データの神経信号化）を施すために、上記のようにして取得した拡散テンソルデータD(l,m,k)を用いて、例えば、拡散テンソルデータD(l,m,k)の固有方程式
【数８】
2007043462.tif 000010
【００４１】
を満たす固有ベクトル
【数９】
2007043462.tif 000011
【００４２】
の中で最大固有値λMに対応する固有ベクトル
【数１０】
2007043462.tif 000012
【００４３】
を算出する（α=1,2,3）。説明を簡単にするために、以降においては、｜λ1｜≧｜λ2｜≧｜λ3｜とし、最大固有値であるλ1に対応する固有ベクトル
【数１１】
2007043462.tif 000013
【００４４】
を
【数１２】
2007043462.tif 000014
【００４５】
と書くことにする。また、該固有ベクトル
【数１３】
2007043462.tif 000015
【００４６】
のノルムは１に規格化されているものとする。【００４７】
次いで、同ステップにおいて、データ前処理手段３は、まず、あるボクセル(l,m,k)に対して、該ボクセル(l,m,k)と隣接する２６個のボクセル（３次元）の中から該ボクセル(l,m,k)上の固有ベクトル
【数１４】
2007043462.tif 000016
【００４８】
との内積の絶対値が最大となる方向ベクトルに対応する隣接ボクセルを求める。この条件を満たす隣接ボクセルを、例えば、(l+1,m+1,k)とする。そして、このようにして求められたボクセルと該ボクセル(l,m,k)を中心として点対称の関係にある隣接ボクセルを求める。この場合は、ボクセル(l-1,m-1,k)である。この状況を図９に示す。上述したように、ボクセル(l+1,m+1,k)及び(l-1,m-1,k)は、ボクセル(l,m,k)上の固有ベクトル
【数１５】
2007043462.tif 000017
【００４９】
との内積が最大となるボクセルであるから、ボクセル(l,m,k)との結合が最も強いと考えられる。なぜならば、図９（ｃ）の場合であれば、ボクセル(l,m,k)を中心としてボクセル(l+1,m+1,k)及び(l-1,m-1,k)は同じ走行方向に走る神経線維上のボクセルと見なすことができるからである。【００５０】
次に、同ステップにおいて、データ前処理手段３は、以下の３通りの状況に応じて、ボクセル(l,m,k)上の脳機能情報の時系列データρ’(l,m,k,i)の書き換えを行う。【００５１】
（状況１）予め取得された所定の脳の解剖学的データ（MRI装置５０により取得されたT1強調画像）から、ボクセル(l+1,m+1,k)及び(l-1,m-1,k)が共に灰白質（神経細胞の細胞体：図中のG）と見なされる場合（図９（Ａ）を参照）：
この場合には、ボクセル(l,m,k)上の脳機能情報の時系列データρ’(l,m,k,i)の値を、
【数１６】
2007043462.tif 000018
【００５２】
と書き換える。【００５３】
（状況２）予め取得された所定の脳の解剖学的データ（MRI装置５０により取得されたT1強調画像）から、ボクセル(l+1,m+1,k)（又は、ボクセル(l-1,m-1,k)）が灰白質と見なされる場合（図９（Ｂ）を参照）：
この場合には、ボクセル(l,m,k)上の脳機能情報の時系列データρ’(l,m,k,i)の値を、
【数１７】
2007043462.tif 000019
【００５４】
と書き換える。【００５５】
（状況３）予め取得された所定の脳の解剖学的データ（MRI装置５０により取得されたT1強調画像）から、ボクセル(l+1,m+1,k)及び(l-1,m-1,k)が共に灰白質とは見なされない（つまり、白質（神経の神経繊維：図中のW）と見なされる）場合（図９（Ｃ）を参照）：
この場合は、ボクセル(l,m,k)上の脳機能情報の時系列データρ’(l,m,k,i)の値を、
【数１８】
2007043462.tif 000020
【００５６】
と書き換える。【００５７】
このような処理を全ボクセルに対して行う。その結果、図９（Ａ），（Ｂ）の下図に示したように、灰白質と見なされるボクセル（灰白質ボクセル）上の脳機能情報の時系列データの値を白質と見なされる隣接ボクセル（白質ボクセル）上の脳機能情報の時系列データの値に反映（伝播）させることができる。また、図９（Ｃ）の下図に示したように、白質と見なされるボクセル同士が隣接している場合には、これらのボクセル上の脳機能データの値が平滑化されることになる。このような前処理を施すことによって、以降のステップの効果を際立たせることができる。【００５８】
なお、このような前処理（ｄ）の仕方は、上記した例に限らず、例えば、灰白質ボクセルに隣接し、且つ該灰白質ボクセルの脳機能情報の時系列データの値を受け取った白質ボクセルが、さらに、隣接する白質ボクセルに対して上記と同様の手順によってその値を順次伝播させるようにしてもよい。また、上記のように灰白質ボクセルの脳機能情報の時系列データの値を白質ボクセルに渡す際に、拡散テンソルデータD(l,m,k)に応じて重み付けをした値を渡すような構成にしてもよい。さらに、上記の手順で書き換えられた白質ボクセルの値と元々の白質ボクセルの値のうちの大きい方を白質ボクセルの値としてもよい。【００５９】
次に、ステップＳ３０において、ボクセル間結合度算出手段４は、上記の拡散テンソルデータ
【数１９】
2007043462.tif 000021
【００６０】
に基づき、隣接ボクセル間の結合度ベクトル
【数２０】
2007043462.tif 000022
【００６１】
を、例えば、
【数２１】
2007043462.tif 000023
【００６２】
から算出する。ここで、
【数２２】
2007043462.tif 000024
【００６３】
の各成分は、例えば、
【数２３】
2007043462.tif 000025
【００６４】
と定義されるのものであり、以降において、これを隣接ボクセル間の平均拡散度ベクトルと称する。αは
【数２４】
2007043462.tif 000026
【００６５】
を満たすパラメータであり、
【数２５】
2007043462.tif 000027
【００６６】
は
【数２６】
2007043462.tif 000028
【００６７】
を満たす定数ベクトルである。【００６８】
式（１３）中の
【数２７】
2007043462.tif 000029
【００６９】
は、プロトンの拡散度を表すベクトル（以降、これを拡散度ベクトルと称する）であり、第１成分Dl(l,m,k)は、位置(l,m,k)のボクセル内のX軸方向へのプロトンの拡散度を表し、第２成分Dm(l,m,k)は、位置(l,m,k)のボクセル内のY方向の拡散度を表し、第３成分Dk(l,m,k)は、位置(l,m,k)のボクセル内のZ方向の拡散度を表す。ここで、X，Y，Z方向とは図７に示した方向と同一の方向を想定している。なお、式（９）においては、拡散テンソルが対称テンソルであることを前提としている。【００７０】
このような意味合いを有する拡散度ベクトル
【数２８】
2007043462.tif 000030
【００７１】
を用いて、式（１３）を介して式（１１）のように構成された隣接ボクセル間の結合度ベクトル
【数２９】
2007043462.tif 000031
【００７２】
の各成分Cl(l,m,k)，Cm(l,m,k)，Ck(l,m,k)はそれぞれ、位置(l,m,k)のボクセルと位置(l+1,m,k)のボクセルとの結合の度合い、位置(l,m,k)のボクセルと位置(l,m+1,k)のボクセルとの結合の度合い、位置(l,m,k)のボクセルと位置(l,m,k+1)のボクセルとの結合の度合いを表す。なお、本実施形態においては、斜め方向に隣接するボクセル間（たとえば、ボクセル(l,m,k)とボクセル(l+1,m+1,k)）の結合度は直接的には考慮していないが、直交する方向のボクセル（この場合、ボクセル(l+1,m,k)やボクセル(l,m+1,k)）を介して間接的に考慮されている。もちろん、斜め方向に隣接するボクセル間の結合度を直接考慮するような構成にしてもよい。【００７３】
さらに、拡散度ベクトル
【数３０】
2007043462.tif 000032
【００７４】
としては、例えば、
（Ｉ）楕円体モデルにおける
【数３１】
2007043462.tif 000033
【００７５】
を採用する。ここで、｜D｜は、式（９）で表される行列の行列式を表す。また、例えば、
（II）式（９）で表される行列の対角成分Dll(l,m,k)，Dmm(l,m,k)，Dkk(l,m,k)を用いて、
【数３２】
2007043462.tif 000034
【００７６】
とすることもできる。また、例えば、
（III）FA（Fractional Anisotropy）：
【数３３】
2007043462.tif 000035
【００７７】
やRA（Relational Anisotropy）：
【数３４】
2007043462.tif 000036
【００７８】
に代表される拡散の異方性を表す指標を用いて拡散度ベクトル
【数３５】
2007043462.tif 000037
【００７９】
を決めることも可能である。ここで、DAVは、
【数３６】
2007043462.tif 000038
【００８０】
である。その他にも、
（IV）予め取得された所定の脳の解剖学的データ（MRI装置５０により取得されたT1強調画像）から脳の白質ボクセルと灰白質ボクセルとを判断し、白質ボクセルでは上記（Ｉ）の拡散度ベクトルを用い、灰白質ボクセルでは、その他の拡散度ベクトル（例えば、上記（II）又は（III））を用いることも可能である。【００８１】
このように、本発明における拡散度ベクトル
【数３７】
2007043462.tif 000039
【００８２】
としては、MRI装置５０から取得可能な拡散テンソルデータD(l,m,k)から神経線維の走行方向が特定可能な量であればいかなるものを用いてもよい。【００８３】
また、パラメータαに関する不等式（１４）の代わりに、例えば、灰白質ボクセルで大きく、白質ボクセルで小さい値といったようにパラメータαの値を灰白質ボクセルと白質ボクセルとで切替えるように設定をしてもよい。【００８４】
また、式（１１）の結合度ベクトル
【数３８】
2007043462.tif 000040
【００８５】
の代わりに、
【数３９】
2007043462.tif 000041
【００８６】
を用いてもよい。【００８７】
さらに、神経繊維の走行方向を強調するために、式（１１）又は式（２２）の結合度ベクトル
【数４０】
2007043462.tif 000042
【００８８】
の各成分cl(l,m,k)，cm(l,m,k)，ck(l,m,k)をN乗（Nは2以上の自然数）したものを採用してもよい。また、これらの成分の値を白質ボクセルでa倍、灰白質ボクセルでb倍したものを用いてもよい。また、これらの成分の中で最も大きいもののみを例えば上記の方法で強調してもよい。さらに、｛cl(l,m,k)＋cm(l,m,k)＋ck(l,m,k)｝/3以上の値を有する成分のみを例えば上記の方法で強調してもよい。【００８９】
さらに、灰白質ボクセルは隣接する６個のボクセル（６方位）とのみ結合させ、白質ボクセルは固有方程式（５）の最大固有値に対応する固有ベクトルとの内積が最大のボクセル及びその点対称のボクセル（２方位）とのみ結合するように結合度ベクトル
【数４１】
2007043462.tif 000043
【００９０】
を構成してもよい。【００９１】
次に、ステップＳ４０において、データ評価量構成手段５は、ステップＳ２０において前処理された脳機能情報の時系列データρ(l,m,k,i)とステップＳ４０において算出された結合度ベクトル
【数４２】
2007043462.tif 000044
【００９２】
とから、図５（Ａ）に例示したようなデータをノンパラメトリック回帰分析するために、例えば、隣接ボクセルを６個とした場合には、データ評価量Q：
【数４３】
2007043462.tif 000045
【００９３】
を構成する。ここで、iは時間、l,m,kはそれぞれボクセルのX方向、Y方向及びZ方向の位置を表す変数であることは上述した通りである。ただし、実際には、図７に示す立体構成されたボクセルの中には脳（神経）を表さないボクセルもあるので、その領域に関するl，m，kは除外するものとする。また、φ(i)は、図５（Ａ）におけるCLASSを表す変数の時系列データであり、この場合、1（=ON）と0（=OFF）のニ値をとる。本実施形態に係るノンパラメトリック回帰分析では、図５（Ａ）に示したような脳機能情報の時系列データρ(l,m,k,i)を説明変数、CLASSを表す変数の時系列データφ(i)を被説明変数とし、これらの変数の間に、
【数４４】
2007043462.tif 000046
【００９４】
を回帰係数とした線形回帰式：
【数４５】
2007043462.tif 000047
【００９５】
を想定する。
【数４６】
2007043462.tif 000048
【００９６】
及び
【数４７】
2007043462.tif 000049
【００９７】
はそれぞれ推定解を表し、κは隣接ボクセル間の連続性を表す重みパラメータである。このように、式（２３）は隣接ボクセル間の連続性に基づいた評価量を構成したが、その代わりに、隣接ボクセル間の滑らかさに基づいた評価量を構成してもよい。【００９８】
本実施形態におけるデータ解析の手法としてノンパラメトリック回帰分析を用いる根拠は、被説明変数φ(i)の数が説明変数ρ(l,m,k,i)の数に比べて極端に少ないことにある。図５（Ａ）に示す例では、説明変数ρ(l,m,k,i)はボクセルの数64×64×64＝262144個あるのに対して、被説明変数φ(i)の測定値の数は18個である。通常でも、被説明変数の測定値の数は数10であり、100を越えることもある。このような状況においては通常の回帰分析ができないので、本実施形態ではノンパラメトリック回帰分析を用いる（非特許文献６：“Spline Smoothing and Nonparametric Regression”, R. L. Eubank, Marcel Dekker, Newyork, 1988）。【００９９】
通常、ノンパラメトリック回帰分析では、被説明変数φ(l,m,k,i)に連続性を仮定することによって、次の評価量Q’：
【数４８】
2007043462.tif 000050
【０１００】
を最小にするような
【数４９】
2007043462.tif 000051
【０１０１】
を求める。ここで、右辺第１項は、各時刻iにおける被説明変数φ(i)の推定解
【数５０】
2007043462.tif 000052
【０１０２】
からの誤差を表す項であり、第２項が、被説明変数φ(i)の時間方向に対する連続性を評価するための項である。しかしながら、図５（Ａ）に示す例からわかるように被被説明変数φ(i)の値に連続性を仮定することはできない。そこで、本実施形態では、このような連続性の制約を被説明変数φ(i)に対してではなく、説明変数ρ(l,m,k,i)に対して課すように拡張したノンパラメトリック回帰分析（非特許文献７：“脳機能画像のノンパラメトリック回帰分析”,月本  洋  他,電子情報通信学会誌D-II，Vol. J84, No.12，pp2623-2633，2001年12月）に基づいて、式（２３）で表されるデータ評価量Qを構成する。【０１０３】
ここで、式（２３）の意味について説明する。右辺第１項は、式（２５）と同じく、各時刻iにおける被説明変数φ(i)の推定解
【数５１】
2007043462.tif 000053
【０１０４】
からの誤差を表す項である。第２項〜第４項は、隣接するボクセルが被説明変数φ(i)に与える影響の連続性を表す項であり、この連続性は隣接する説明変数ρ(l,m,k,i)が被説明変数に与える影響、すなわち、式（２４）の回帰係数
【数５２】
2007043462.tif 000054
【０１０５】
の連続性、で評価できる。【０１０６】
さらに、本実施形態においては、第２項〜第４項において、ボクセル間結合度算出手段４で算出された結合度ベクトル
【数５３】
2007043462.tif 000055
【０１０７】
の各成分Cl(l,m,k)，Cm(l,m,k)，Ck(l,m,k)をそれぞれX方向、Y方向、Z方向に対する隣接ボクセル間の連続性に対する重み付け因子として導入している。式（１１）に示したように、これらの成分はそれぞれ、平均拡散度ベクトル
【数５４】
2007043462.tif 000056
【０１０８】
の各成分D’l(l,m,k)，D’m(l,m,k)，D’k(l,m,k)の関数となっているので、結果として、式（２３）のデータ評価量Qは神経線維（神経細胞）の走行方向（プロトンの拡散の異方性）を反映していることになる。【０１０９】
次に、ステップＳ５０において、データ解析手段６は、ステップＳ４０において構成されたデータ評価量Qの極値を求めるための計算を行う。本実施形態においては、式（２３）のデータ評価量Qを最小にするような回帰係数
【数５５】
2007043462.tif 000057
【０１１０】
を算出する。ここでは、簡単のため、式（２３）を
【数５６】
2007043462.tif 000058
【０１１１】
のように書き換えたデータ評価量Q”について説明する。ここで、nは立体構成された全ボクセルに対して所定の番号付けを施すことによって１次元配列し直したボクセルの位置を表す正の離散変数であり、３変数表示(l,m,k)と１対１対応している。さらに、簡単のため、ここではL=M=Kとしたので、全ボクセル数はL3となる。また、C(n)は、このように１次元に配列し直した際に、ボクセルnと該ボクセルnと隣接するボクセルn’との間の結合度を表す係数であり、この係数は、式（１１）の結合ベクトル
【数５７】
2007043462.tif 000059
【０１１２】
から一意に定まる。例えば、図３（Ｂ）に示した２次元スライス画像Skに話を限定した場合、n=130のボクセルと隣接するボクセルはn=127,129,131,257の４つとなる。ステップＳ３０で算出された結合度ベクトル
【数５８】
2007043462.tif 000060
【０１１３】
により、ボクセルn=127とボクセルn=130との結合度、ボクセルn=129とボクセルn=130との結合度、ボクセルn=131とボクセルn=130との結合度、ボクセルn=257とボクセルn=130との結合度はそれぞれ、例えば、2，2，5，5のように決まる。つまり、この場合、
【数５９】
2007043462.tif 000061
【０１１４】
となる。【０１１５】
このとき、回帰係数
【数６０】
2007043462.tif 000062
【０１１６】
は、最小２乗法を用いて、
【数６１】
2007043462.tif 000063
【０１１７】
のように求めることができる。ここで、Xはρ(n,i)を成分とするL3×I行列Xと、
【数６２】
2007043462.tif 000064
【０１１８】
はφ(i)を成分とするI次元ベクトル、
【数６３】
2007043462.tif 000065
【０１１９】
は
【数６４】
2007043462.tif 000066
【０１２０】
を成分とするL3次元ベクトル、κOは交差検証法（Cross Validation）で求めたκの最適値、すなわち、
【数６５】
2007043462.tif 000067
【０１２１】
を最小にするようなκである。ここで、
【数６６】
2007043462.tif 000068
【０１２２】
はi番目の事例を除いて回帰係数を計算し、その回帰係数を用いてi番目の事例の測定値を予測した値である。【０１２３】
次に、ステップＳ６０において、画像生成手段７は、上記のようにして算出された回帰係数
【数６７】
2007043462.tif 000069
【０１２４】
を、例えば、回帰係数の値が0（すなわち、無相関）のボクセルをグレイスケールの基準として、回帰係数の値が正（すなわち、正の相関）のボクセルを白方向の色調とし、回帰係数の値が負（すなわち、負の相関）のボクセルを黒方向の色調とするような色調を施した画像を生成する。その他にも、グレイスケール上で取り得るボクセルの値の頻度のヒストグラムにおいて、上位５％の値以上を表示させるようにしてもよいし、所定の規則に従ってフルカラー画像を生成してもよい。【０１２５】
そして、ステップＳ７０において、画像表示手段８は、上記のようにして生成された画像を立体表示する。【０１２６】
このように、第１の実施形態においては、MRI装置から取得した脳機能情報の時系列データに基づく誤差の評価項に加えて、同じくMRI装置から取得した拡散テンソルデータに基づく空間方向の連続性の評価項を加えた評価量を構成した。そして、この評価項をノンパラメトリック線形回帰分析し、その結果得られた回帰係数を画像表示の対象としたことによって、脳の活性部位間の結合構造をも考慮した脳の活性部位の特定が可能となる。【０１２７】
なお、本実施形態においては、ステップＳ２０において通常の前処理（ａ）〜（ｃ）以外に前処理（ｄ）のプロセスを設けたが、前処理（ｄ）をはずしてもよい。このことは、以降の実施形態においても同様である。【０１２８】
なお、第１の実施形態では、データ評価量として式（２３）を構成したわけだが、本発明の目的を達成するにあたり、データ解析の手法に応じて、その他にも各種の評価量を構成することができる。そこで、次に、上記した第１の実施形態の幾つかの変形例について説明するが、重複を避けるため、以降においては第１の実施形態と異なる部分についてのみ説明を加える。【０１２９】
＜第１の実施形態の変形例１＞
第１の実施形態では、ステップＳ５０において、データ解析手段６は、式（１９）の線形回帰式を用いることによりノンパラメトリック回帰分析をしたわけだが、より厳密には、回帰式を２次関数、３次関数等の一般のｎ次関数で行うことも考えられる。このような非線形回帰分析を行う際には、ニューラルネットワーク（Neural Network）モデルを用いることができる。典型的なニューラルネットワークモデルとして、入力層、中間層、出力層から成る階層型ニューラルネットワークモデルがある。【０１３０】
本変形例では、ステップＳ４０において、データ評価量構成手段５は、ニューラルネットワークモデルにおける通常の学習の評価量である２乗誤差の項に対して、ボクセル間結合度算出手段４で算出された結合度ベクトル
【数６８】
2007043462.tif 000070
【０１３１】
に基づく隣接ボクセル間の連続性を評価する評価量を加えた評価量を構成する。そして、ステップＳ５０において、データ解析手段６は、このようにして構成された評価量の最小値を求めるために、誤差逆伝播法（Error Back Propagation Method）等を用いて解析した後に、ステップＳ６０において、画像生成手段７は、入力（脳機能情報の時系列データ）と出力の間の感度（入力を変化させたときに出力に与える影響）を表示画像として第１の実施形態と同様にして生成する。そして、ステップＳ７０において、画像表示手段８は、上記のようにして生成された画像を立体表示する。【０１３２】
これによって、第１の実施形態と同様の効果が得られる。【０１３３】
＜第１の実施形態の変形例２＞
第１の実施形態では、ステップＳ４０において、データ評価量構成手段５は、脳機能の時系列データρ(l,m,k,i)を説明変数、φ(i)を被説明変数として式（２３）の評価量Qを構成したわけだが、本変形例では、その逆、つまり、φ(i)を説明変数、ρ(l,m,k,i)を被説明変数として、ボクセル(l,m,k)毎に、例えば、
【数６９】
2007043462.tif 000071
【０１３４】
といった評価量Q”’を構成する。ここで、
【数７０】
2007043462.tif 000072
【０１３５】
である。【０１３６】
そして、ステップＳ５０において、データ解析手段６は、このようにして構成された評価量Q”’の最小値を求めるために、最適なκを交差検証法を用いて決定し、GA（Genetic Algorithm:遺伝的アルゴリズム）等の確率的探索アルゴリズムによって回帰係数
【数７１】
2007043462.tif 000073
【０１３７】
を求める。【０１３８】
次いで、ステップＳ６０において、画像生成手段７は、上記の解析結果を各種検定した際に有意であると判断されたボクセルの脳機能データρ(l,m,k)に基づく値を表示画像として生成する。そして、ステップＳ７０において、画像表示手段８は、上記のようにして生成された画像を立体表示する。【０１３９】
これによって、第１の実施形態と同様の効果が得られる。【０１４０】
＜第１の実施形態の変形例３＞
第１の実施形態の変形例１では、ステップＳ４０において、データ評価量構成手段５は、脳機能の時系列データρ(l,m,k,i)を説明変数、φ(i)を被説明変数として式（２３）の評価量を構成したわけだが、本変形例では、その逆、つまり、φ(i)を説明変数、ρ(l,m,k,i)を被説明変数として、階層型ニューラルネットワークモデルを用いて、非線形回帰分析を行うこともできる。ただし、この場合には、ステップＳ５０において極値を求める際に誤差逆伝播法が使えないので、その代わりに、データ解析手段６は、最適なκを交差検証法を用いて決定し、GA等の確率的探索アルゴリズムによってデータ解析を行う。そして、ステップＳ６０において、画像生成手段７は、上記の解析結果を各種検定した際に有意であると判断されたボクセルの脳機能データρ(l,m,k)に基づく値を表示画像として生成する。そして、ステップＳ７０において、画像表示手段８は、上記のようにして生成された画像を立体表示する。【０１４１】
これによって、第１の実施形態と同様の効果が得られる。【０１４２】
＜第１の実施形態の変形例４＞
第１の実施形態では、ステップＳ５０において、データ解析手段６は、式（２４）の線形回帰式を用いることによりノンパラメトリック回帰分析をしたわけだが、データ解析の手法として、判別分析（非特許文献８：“多変量解析のはなし”，有馬哲，石村貞夫， 東京図書，1987）、数量化II類（非特許文献９）、決定木（非特許文献１０：“AIによるデータ解析”，J.R.キンラン，トッパン，1995）、サポートベクターマシン（非特許文献１１：“サポートベクターマシン入門”，Nello Cristianini and John Shawe‐Taylor，共立出版，2005）等の分類的解析手法を用いることもできる。【０１４３】
いずれの場合も、例えば、説明変数を脳機能情報の時系列データρ(l,m,k,i)、外的基準をタスクとレストの時系列データとして、各種の判別関数、決定木等を用いてデータ分類する。本変形例では、その際に、ステップＳ４０において、データ評価量構成手段５は、各解析手法における通常の評価量に対して、ステップＳ３０で算出した結合度ベクトル
【数７２】
2007043462.tif 000074
【０１４４】
を組み込んだ評価量を構成する。そして、ステップＳ５０において、データ解析手段６は、各解析手法を用いてデータ解析を行う。上記のように評価量を変更すると、通常の解法が使えなくなる場合があるので、その際には、データ解析手段６は、最適なκを交差検証法を用いて決定し、GA等の確率的探索アルゴリズムによってデータ解析を行えばよい。【０１４５】
そして、ステップＳ６０において、画像生成手段７は、ステップＳ５０で得られた分類データを各種検定した際に有意であると判断されたボクセルの脳機能データρ(l,m,k)に基づく値を表示画像として生成し、ステップＳ７０において、画像表示手段８は、上記のようにして生成された画像を立体表示する。【０１４６】
これによって、第１の実施形態と同様の効果が得られる。【０１４７】
＜第１の実施形態の変形例５＞
第１の実施形態では、ステップＳ５０において、データ解析手段６は、式（２４）の線形回帰式を用いることによりノンパラメトリック回帰分析をしたわけだが、データ解析の手法として、ICA（Independent Component Analysis：独立成分分析）を用いることもできる（非特許文献１２：“詳解独立成分分析 信号解析の新しい世界”,  Aapo Hyvarinen et al.，東京電機大学出版局  2005）。【０１４８】
本変形例では、ステップＳ４０において、データ評価量構成手段５は、ボクセル間結合度算出手段４で算出された結合度ベクトル
【数７３】
2007043462.tif 000075
【０１４９】
に基づき評価量を構成する。独立成分分析による脳機能情報の時系列データρ(l,m,k,i)の解析での「空間的な独立」という条件を、本変形例では「拡散テンソル情報にもとづく結合度を考慮した空間的な独立」というように変更する。独立成分分析の具体的手法はいくつかあるが、例えば、相互相関を最小化する手法の場合には、本変形例では、相互相関を最小化するのではなく、結合度に基づく評価量に近づけるというようにする。これにより、結合度の大きいボクセル間の相互相関は結合度に基づく評価量に近いような「独立度」になり、結合度の小さい（もしくは、結合度のない）ボクセル間は相互相関を最小にするというようになる。より具体的には、例えば、相関行列の非対角成分の２乗和を最小化するような評価関数中の２乗和の各項に、結合度に基づく評価量による重み係数をつける。この重み係数は、結合度が大きいほど小さくなるように構成する。【０１５０】
次に、ステップＳ５０において、データ解析手段６は、上記のようにして構成された評価量を考慮した独立成分分析を行う。その際、最適なκを交差検証法を用いて決定し、GA等の確率的探索アルゴリズムによってデータ解析を行う。【０１５１】
次に、ステップＳ６０において、画像生成手段７は、上記のようにして得られた独立成分を表示画像として生成する。ただし、上記の解析の結果得られる独立成分の数は多いので、全行ベクトルの中からタスク／レストと最も相関の高い行ベクトルのみを抽出する処理を行う。【０１５２】
そして、ステップＳ７０において、画像表示手段８は、上記のようにして生成された画像を立体表示する。【０１５３】
これによって、第１の実施形態と同様の効果が得られる。【０１５４】
なお、上記した独立成分分析以外にも、脳機能情報の時系列データから所定の主たる成分を抽出する手法として、主成分分析、因子分析、数量化III類等があるが、これらの解析手法についても、基本的には、上記と同様にして解析することができる。【０１５５】
[第２の実施形態]
図１０は、本発明の第２の実施形態に係る脳機能解析装置の概略構成を示したブロック図である。第２の実施形態の脳機能解析装置２０は、第１の実施形態の脳機能解析装置１０におけるボクセル間結合度算出手段４、データ評価量構成手段５、及びデータ解析手段６の代わりに、ボクセル間結合度算出手段２４、データ評価量構成手段２５、データ平滑化手段２００、及びデータ解析手段２６を設けた構成を成している。第１の実施形態と同一の構成には同一の符号を付し、以降においては、繰り返しの説明を省略する。なお、本実施形態においては、脳機能解析装置２０を、本実施形態の主要構成手段２０Ａと画像表示手段８及び記憶手段９とを一体型としたディスプレイコンソール形式としたが、画像表示手段８、或いは、画像表示手段８及び記憶手段９をそれぞれ独立した画像表示装置及び記憶装置として主要構成手段２０Ａから切り離した構成としてもよい。【０１５６】
ボクセル間結合度算出手段２４は、拡散テンソルデータ取得手段２で取得された拡散テンソルデータD(l,m,k)（式（９））から第１の実施形態に詳述したようにして平均拡散度ベクトル
【数７４】
2007043462.tif 000076
【０１５７】
（式（１２））を算出した後、該平均拡散度ベクトル
【数７５】
2007043462.tif 000077
【０１５８】
の各成分D’l(l,m,k)，D’m(l,m,k)，D’k(l,m,k)の大きさに応じてデータ前処理手段３で前処理された脳機能情報の時系列データρ(l,m,k,i)を平滑化する際の隣接ボクセル間の平滑度を決める結合度ベクトル
【数７６】
2007043462.tif 000078
【０１５９】
を算出する。【０１６０】
データ平滑化手段２００は、ボクセル間結合度算出手段２４で算出された結合度ベクトル
【数７７】
2007043462.tif 000079
【０１６１】
を用いて、データ前処理手段３で前処理された脳機能情報の時系列データρ(l,m,k,i)を平滑化するための加重平均フィルターを構成し、該加重平均フィルターを用いて、データ前処理手段３で前処理された脳機能情報の時系列データρ(l,m,k,i)を平滑化する。【０１６２】
データ解析手段２６は、データ平滑化手段２００で平滑化されたデータをSPM等の手法を用いてデータ解析を行う。【０１６３】
図１１は、図１０に示した第２の実施形態に係る脳機能解析装置２０によって行われるデータ解析の処理手順の一例を示したフローチャートである。【０１６４】
ステップＳ１００において、脳機能データ取得手段１は、MRI装置５０で計測された脳機能情報の原時系列データρ’(l,m,k,i)を取得する。同ステップにおいて、拡散テンソル情報取得手段２は、同じくMRI装置５０で計測された拡散テンソルデータD(l,m,k)を取得する。【０１６５】
次に、ステップＳ１１０において、データ前処理手段３は、ステップＳ１００において取得した脳機能情報の原時系列データρ’(l,m,k,i)に対して、第１の実施形態で詳述した前処理（ａ）〜（ｄ）を施すことによって、脳機能情報の時系列データρ(l,m,k,i)を生成する。【０１６６】
次に、ステップＳ１２０において、ボクセル間結合度算出手段２４は、まず、ステップＳ１００において取得した拡散テンソルデータD(l,m,k)から第１の実施形態で詳述したようにして平均拡散度ベクトル
【数７８】
2007043462.tif 000080
【０１６７】
を算出する。次いで、該平均拡散度ベクトル
【数７９】
2007043462.tif 000081
【０１６８】
の各成分D’l(l,m,k)，D’m(l,m,k)，D’k(l,m,k)の大きさに応じてデータ前処理手段３で前処理された脳機能情報の時系列データρ(l,m,k,i)を平滑化する際の隣接ボクセル間の平滑度を決める結合度ベクトル 
【数８０】
2007043462.tif 000082
【０１６９】
を、例えば、
【数８１】
2007043462.tif 000083
【０１７０】
として算出する。【０１７１】
次に、ステップＳ１３０において、データ平滑化手段２６は、ステップＳ１２０において算出された結合度ベクトル
【数８２】
2007043462.tif 000084
【０１７２】
を用いて、ステップＳ１１０において前処理された脳機能情報の時系列データρ(l,m,k,i)を、例えば、
【数８３】
2007043462.tif 000085
【０１７３】
のように異方的に平滑化する。ここで、w及びWは重み因子であり、
【数８４】
2007043462.tif 000086
【０１７４】
を満たすように規格化する。【０１７５】
図１２は、本実施形態に係るデータ平滑化の手法を例示したものであり、ある２次元スライス画像Skの位置(l,m)を中心とした９つのボクセル上の脳機能情報の時系列データρk(l,m,i)を表している。この例では、ρk(l,m,i)=5，ρk(l-1,m,i)=ρk(l-1,m-1,i)=ρk(l,m-1,i)=ρk(l,m+1,i)=ρk(l+1,m,i)=2，ρk(l+1,m-1,i)=ρk(l-1,m+1,i)=ρk(l+1,m+1,i)=3であり、このうち、ρk(l,m,i)の平滑化に用いられるのは、ρk(l,m-1,i)，ρk(l,m+1,i)，ρk(l-1,m,i)，ρk(l+1,m,i)の４つである。さらに、D’k(l,m,k)=0（つまり、Z方向へのプロトンの拡散がないということ）、且つ、D’l(l,m,k)=0.8，D’m(l,m,k)=0.2，w=0.4，W=0.6を満たすような状況を想定する。このとき、結合度ベクトル
【数８５】
2007043462.tif 000087
【０１７６】
は、式（３１）から
【数８６】
2007043462.tif 000088
【０１７７】
となり、位置(l,m,k)のボクセル上の脳機能情報の時系列データρ(l,m,k,i)=5は、式（３６）から、
【数８７】
2007043462.tif 000089
【０１７８】
のように平滑化されることがわかる。この処理を該２次元スライス画像Skの全ボクセル上の脳機能情報の時系列データρ(l,m,k,i)に施すことによって、神経繊維の走行方向に沿った方向（つまり、平均拡散度ベクトル
【数８８】
2007043462.tif 000090
【０１７９】
の成分D’l(l,m,k)，D’m(l,m,k)，D’k(l,m,k)のうち、最も大きい成分の方向：この例では、X方向）からの寄与が大きいように平滑化することができる。【０１８０】
また、これ以外にも、式（５）の箇所で説明したように、最大固有値λMに対応する固有ベクトル
【数８９】
2007043462.tif 000091
【０１８１】
の向きと合致する２個のボクセルとのみ平滑化してもよい。【０１８２】
次に、ステップＳ１４０において、データ解析手段２６は、ステップＳ１４０において平滑化された脳機能情報の時系列データ
【数９０】
2007043462.tif 000092
【０１８３】
に対して、SPM等の解析手法を用いてデータ解析を施す。【０１８４】
次に、ステップＳ１５０において、画像生成手段７は、データ解析手段２６によって解析された結果をグレイスケールやフルカラー画像等の色調画像を生成する。【０１８５】
そして、ステップＳ１６０において、画像表示手段８は、上記のようにして生成された画像を立体表示する。【０１８６】
このように、第２の実施形態においては、MRI装置から取得した脳機能情報の時系列データを、同じくMRI装置から取得した拡散テンソルデータから算出された結合度ベクトルに基づいて平滑化処理したことによって、脳の活性部位間の結合構造をも考慮した脳の活性部位の特定が可能となる。【０１８７】
なお、式（３１）の結合度ベクトル
【数９１】
2007043462.tif 000093
【０１８８】
や式（３２）の平滑化処理のアルゴリズムはあくまでも一例であって、その他の実施例も可能である。さらに、本実施形態におけるデータ平滑化処理の手法は、特定のデータ解析手法に依存しないことは上記の説明から明らかである。【０１８９】
＜第２の実施形態の変形例＞
第２の実施形態においては、ステップＳ１３０において脳機能情報の時系列データρ(l,m,k,i)に対して拡散テンソルデータD(l,m,k)から算出された結合度ベクトル
【数９２】
2007043462.tif 000094
【０１９０】
を用いて平滑化処理を施した後に、ステップＳ１４０においてSPM等の解析手法を用いてデータ解析を行ったが、本発明の目的を達成するにあたって、これらのステップを逆にするような手法も考えられる。そこで、次に、上記した第２の実施形態の変形例について説明する。なお、重複を避けるため、以降においては第２の実施形態と異なる部分についてのみ説明を加える。【０１９１】
図１３は、第２の実施形態の変形例に係る脳機能解析装置の概略構成を示したブロック図である。本変形例に係る脳機能解析装置２０’の構成は、第２の実施形態の構成と基本的には同じであるが、データ平滑化手段２００とデータ解析手段２６とが入れ替わった構成になっている。これによって、脳機能解析装置２０’によって行われるデータ解析の処理手順が第２の実施形態と異なる。なお、本実施形態においては、脳機能解析装置２０’を、本実施形態の主要構成手段２０Ａ’と画像表示手段８及び記憶手段９とを一体型としたディスプレイコンソール形式としたが、画像表示手段８、或いは、画像表示手段８及び記憶手段９をそれぞれ独立した画像表示装置及び記憶装置として主要構成手段２０Ａ’から切り離した構成としてもよい。【０１９２】
図１４は、図１３に示した脳機能解析装置２０’によって行われるデータ解析の処理手順の一例を示したフローチャートである。本変形例では、ステップＳ１２１において、まず、データ解析手段２６’は、データ前処理手段３によって前処理された脳機能情報の時系列データρ(l,m,k,i)をSPM等の解析手法を用いてデータ解析する。次に、ステップＳ１３１において、ボクセル間結合度算出手段２４は、第２の実施形態と同様にして、式（３１）で示したような結合度ベクトル
【数９３】
2007043462.tif 000095
【０１９３】
を算出する。次いで、ステップＳ１４１において、データ平滑化手段２００’は、データ解析手段２６’によって解析された結果である脳機能データρ(l,m,k)に対して、第２の実施形態と同様にして、式（３２）を用いて平滑化処理を行う。次いで、ステップＳ１５０において、画像生成手段７は、上記のようにして得られた平滑化処理後のデータ
【数９４】
2007043462.tif 000096
【０１９４】
に対して、所定の処理規則に従って、グレイスケールやフルカラー画像等の色調画像を生成する。そして、ステップＳ１６０において、画像表示手段８は、上記のようにして生成された画像を立体表示する。【０１９５】
このようにして、本変形例においても、第２の実施形態と同様の効果が得られる。【０１９６】
[第３の実施形態]
図１５は、本発明の第３の実施形態に係る脳機能解析装置の概略構成を示したブロック図である。第３の実施形態の脳機能解析装置３０は、第１の実施形態の脳機能解析装置１０におけるボクセル間結合度算出手段４、データ評価量構成手段５、及びデータ解析手段６の代わりに、ボクセル間結合度算出手段３４、クラスタリング手段３００、及びデータ解析手段３６を設けた構成を成している。第１の実施形態と同一の構成には同一の符号を付し、以降においては、繰り返しの説明を省略する。なお、本実施形態においては、脳機能解析装置３０を、本実施形態の主要構成手段３０Ａと画像表示手段８及び記憶手段９とを一体型としたディスプレイコンソール形式としたが、画像表示手段８、或いは、画像表示手段８及び記憶手段９をそれぞれ独立した画像表示装置及び記憶装置として主要構成手段３０Ａから切り離した構成としてもよい。【０１９７】
ボクセル間結合度算出手段３４は、拡散テンソルデータ取得手段２で取得された拡散テンソルデータD(l,m,k)（式（９））から第１の実施形態で詳述したようにして平均拡散度ベクトル
【数９５】
2007043462.tif 000097
【０１９８】
（式（１２））を算出した後、該平均拡散度ベクトル
【数９６】
2007043462.tif 000098
【０１９９】
の各成分D’l(l,m,k)，D’m(l,m,k)，D’k(l,m,k)の大きさに応じてボクセル同士をクラスタリングするための指標として結合度ベクトル
【数９７】
2007043462.tif 000099
【０２００】
を算出する。【０２０１】
クラスタリング手段３００は、ボクセル間結合度算出手段３４によって算出された結合度ベクトル
【数９８】
2007043462.tif 000100
【０２０２】
を用いてボクセルをクラスタリングする。【０２０３】
データ解析手段３６は、クラスタリング手段３００でクラスタリングされたボクセル集団（クラスター）を単位として、データ前処理手段３で前処理された脳機能情報の時系列データρ(l,m,k,i)をSPM等の解析手法を用いてデータ解析を行う。【０２０４】
図１６は、図１５に示した第３の実施形態に係る脳機能解析装置３０によって行われるデータ解析の処理手順の一例を示したフローチャートである。【０２０５】
ステップＳ２００において、脳機能データ取得手段１は、MRI装置５０で計測された脳機能情報の原時系列データρ’(l,m,k,i)を取得する。同ステップにおいて、拡散テンソル情報取得手段２は、同じくMRI装置５０で計測された拡散テンソルデータD(l,m,k)を取得する。【０２０６】
次に、ステップＳ２１０において、データ前処理手段３は、ステップＳ２００において取得した脳機能情報の原時系列データρ’(l,m,k,i)に対して、第１の実施形態で詳述した前処理（ａ）〜（ｄ）を施すことによって、脳機能情報の時系列データρ(l,m,k,i)を生成する。【０２０７】
次に、ステップＳ２２０において、ボクセル間結合度算出手段３４は、まず、ステップＳ２００において取得した拡散テンソルデータD(l,m,k)から第１の実施形態で詳述したようにして平均拡散度ベクトル
【数９９】
2007043462.tif 000101
【０２０８】
を算出する。次いで、該平均拡散度ベクトル
【数１００】
2007043462.tif 000102
【０２０９】
の各成分D’l(l,m,k)，D’m(l,m,k)，D’k(l,m,k)の大きさに応じてボクセル同士をクラスタリングする際の指標となる結合度ベクトル
【数１０１】
2007043462.tif 000103
【０２１０】
を、例えば、
【数１０２】
2007043462.tif 000104
【０２１１】
として算出する。ここで、
【数１０３】
2007043462.tif 000105
【０２１２】
はクラスタリングの閾値を表す定数ベクトルであり、
【数１０４】
2007043462.tif 000106
【０２１３】
と書ける。【０２１４】
次に、ステップＳ２３０において、クラスタリング手段３００は、ステップＳ２２０において算出された結合度ベクトル
【数１０５】
2007043462.tif 000107
【０２１５】
の各成分Sl(l,m,k)，Sm(l,m,k),Sk(l,m,k)の値が正になるボクセルを、それぞれX方向、Y方向、Z方向に対してクラスタリングする。図１７は、図１５に示したクラスタリング手段３００によって行われるクラスタリングの手法を例示した図であり、ある２次元スライス画像Skの１６個のボクセルを表している。図１７（Ａ）における隣接ボクセル間の矢印の値は、X方向又はY方向の平均拡散度ベクトルの値を表している。式（３７）の定数ベクトル
【数１０６】
2007043462.tif 000108
【０２１６】
の各成分の値をそれぞれ、例えば、hl=hm=3，hk=0とすると、結合度ベクトル
【数１０７】
2007043462.tif 000109
【０２１７】
の成分Sl(l,m,k)，Sm(l,m,k)の各値は、図１７（Ｂ）のようになる。このとき、X方向、Y方向の各方向において、結合度ベクトル
【数１０８】
2007043462.tif 000110
【０２１８】
の成分Sl(l,m,k)，Sm(l,m,k)が共に正となる隣接ボクセルを順次結んでいくことによって、最終的に、点線で囲んだボクセルが一つのクラスターとして形成されることになる。結合度ベクトル
【数１０９】
2007043462.tif 000111
【０２１９】
が平均拡散度ベクトルD’(l,m,k)の関数となっていることから、このようなクラスタリングは神経線維の走行方向を反映したものとなっている。【０２２０】
次に、ステップＳ２４０において、データ解析手段３６は、ステップＳ２３０においてクラスタリングされたボクセル集団（クラスター）毎に、ステップＳ２１０において前処理された脳機能情報の時系列データρ(l,m,k,i)を、SPM等の解析方法を用いてデータ解析する。その際には、例えば、一つのクラスター内での脳機能情報の時系列データを平均化する等の処理が必要となる。【０２２１】
次に、ステップＳ２５０において、画像生成手段７は、ステップＳ２４０においてデータ解析された後のデータ
【数１１０】
2007043462.tif 000112
【０２２２】
に対して、所定の処理規則に従って、グレイスケールやフルカラー画像等の色調画像を生成する。【０２２３】
そして、ステップＳ２６０において、画像表示手段８は、ステップＳ２５０において生成された画像を立体表示する。【０２２４】
このように、第３の実施形態においては、MRI装置から取得した拡散テンソルデータに基づいてボクセルをクラスタリングした後に、同じくMRI装置から取得した脳機能情報の時系列データをクラスター毎に解析したことによって、脳の活性部位間の結合構造をも考慮した脳の活性部位の特定が可能となる。【０２２５】
なお、式（３６）の結合度ベクトル
【数１１１】
2007043462.tif 000113
【０２２６】
はあくまでも一例であってその他の実施例も可能である。また、式（５）の箇所で説明したように、最大固有値λMに対応する固有ベクトル
【数１１２】
2007043462.tif 000114
【０２２７】
の向きと合致する２個のボクセルとのみクラスタリングしてもよい。さらに、本実施形態におけるクラスタリング処理の手法は、特定のデータ解析手法に依存しないことは上記の説明から明らかである。【０２２８】
＜第３の実施形態の変形例１＞
第３の実施形態においては、ステップＳ２３０において、拡散テンソルデータD(l,m,k)から算出された結合度ベクトル
【数１１３】
2007043462.tif 000115
【０２２９】
を用いてボクセルをクラスタリングした後に、ステップＳ２４０においてクラスター毎に脳機能情報の時系列データρ(l,m,k,i)をSPM等の解析手法を用いてデータ解析したが、本発明の目的を達成するにあたって、これらのステップを逆にするような手法も考えられる。そこで、次に、上記した第３の実施形態の変形例について説明する。なお、重複を避けるため、以降においては第３の実施形態と異なる部分についてのみ説明を加える。【０２３０】
図１８は、第３の実施形態の変形例に係る脳機能解析装置の概略構成を示したブロック図である。本変形例に係る脳機能解析装置３０’の構成は、第３の実施形態の構成と基本的には同じであるが、クラスタリング手段３００とデータ解析手段３６とが入れ替わった構成になっている。これによって、脳機能解析装置３０’によって行われるデータ解析の処理手順が第３の実施形態と異なる。なお、本実施形態においては、脳機能解析装置３０’を、本実施形態の主要構成手段３０Ａ’と画像表示手段８及び記憶手段９とを一体型としたディスプレイコンソール形式としたが、画像表示手段８、或いは、画像表示手段８及び記憶手段９をそれぞれ独立した画像表示装置及び記憶装置として主要構成手段３０Ａ’から切り離した構成としてもよい。【０２３１】
図１９は、図１８に示した脳機能解析装置３０’によって行われるデータ解析の処理手順の一例を示したフローチャートである。本変形例では、ステップＳ２２１において、まず、データ解析手段３６’は、データ前処理手段３によって前処理された脳機能情報の時系列データρ(l,m,k,i)をSPM等の解析手法でデータ解析する。次に、ステップＳ２３１において、ボクセル間結合度算出手段３４は、第３の実施形態と同様にして、式（３６）で示したような結合度ベクトル
【数１１４】
2007043462.tif 000116
【０２３２】
を算出する。次いで、ステップＳ２４１において、クラスタリング手段３００’は、データ解析手段３６’によって解析された脳機能データρ(l,m,k)に対して、第３の実施形態と同様にして、図１７に示したようなボクセルのクラスタリングを行う。次いで、ステップＳ２５０において、画像生成手段７は、上記のようにしてクラスタリングされた脳機能データ
【数１１５】
2007043462.tif 000117
【０２３３】
に対して、所定の処理規則に従って、グレイスケールやフルカラー画像等の色調画像を生成する。そして、ステップＳ２６０において、画像表示手段８は、上記のようにして生成された画像を立体表示する。【０２３４】
このようにして、本変形例においても、第３の実施形態と同様の効果が得られる。【０２３５】
＜第３の実施形態の変形例２＞
第３の実施形態の変形例１においては、データ前処理手段３によって前処理された脳機能情報の時系列データρ(l,m,k,i)をSPM等の解析手法でデータ解析した後に、ボクセル間結合度算出手段３４で算出された結合度ベクトル
【数１１６】
2007043462.tif 000118
【０２３６】
を用いてデータ解析後の脳機能データρ(l,m,k)をクラスタリングしたが、データ解析の手法として分類的解析手法を用いることも考えられる。そこで、本変形例では、その一例として、データ解析の手法として分類的解析手法の一つである決定木を用いた場合について説明する。例えば、図４に示したように、説明変数（属性）を脳機能情報の時系列データρ(l,m,k)、外的基準（クラス）をタスクとレストの時系列データとして決定木の手法を用いた場合、不要な属性（ボクセル）は除去され、必要な属性（ボクセル）のみが残る。このようにして残ったボクセルを結合度ベクトル
【数１１７】
2007043462.tif 000119
【０２３７】
を用いてクラスタリングする。例えば、（１）あるボクセルを選び、（２）そのボクセルの周囲で結合度が一定値以上のボクセルをクラスタリングし、（３）（２）の処理を、周囲に結合度が一定値以上のボクセルがなくなるまで続け、（４）ボクセルが残っていれば、あるボクセルを選び、（２）の処理にもどる。【０２３８】
次いで、ステップＳ２５０において、画像生成手段７は、上記のようにしてクラスタリングされた脳機能データ
【数１１８】
2007043462.tif 000120
【０２３９】
に対して、所定の処理規則に従って、グレイスケールやフルカラー画像等の色調画像を生成する。そして、ステップＳ２６０において、画像表示手段８は、上記のようにして生成された画像を立体表示する。【０２４０】
このようにして、本変形例においても、第３の実施形態と同様の効果が得られる。【０２４１】
[第４の実施形態]
図２０は、第４の実施形態に係る脳機能解析装置の概略構成を示したブロック図である。第４の実施形態の脳機能解析装置４０は、第１の実施形態の脳機能解析装置１０におけるボクセル間結合度算出手段４、データ評価量構成手段５、及びデータ解析手段６の代わりに、ボクセル間結合度算出手段４４、データ検定手段４００を設けた構成を成している。第１の実施形態と同一の構成には同一の符号を付し、以降においては、繰り返しの説明を省略する。なお、本実施形態においては、脳機能解析装置４０を、本実施形態の主要構成手段４０Ａと画像表示手段８及び記憶手段９とを一体型としたディスプレイコンソール形式としたが、画像表示手段８、或いは、画像表示手段８及び記憶手段９をそれぞれ独立した画像表示装置及び記憶装置として主要構成手段４０Ａから切り離した構成としてもよい。【０２４２】
ボクセル間結合度算出手段４４は、データ前処理手段３で前処理された脳機能情報の時系列データρ(l,m,k,i)をボクセル毎に各種データ検定（例えば、ｔ検定、ｆ検定、又はｚ検定）する際の基準値（検定値や有意水準など）を調整するための結合度ベクトル
【数１１９】
2007043462.tif 000121
【０２４３】
を、式（１２）の平均拡散度ベクトル
【数１２０】
2007043462.tif 000122
【０２４４】
の関数として、拡散テンソルデータ取得手段２で取得された拡散テンソルデータD(l,m,k)（式（９））から算出する。【０２４５】
データ検定手段４００は、ボクセル間結合度算出手段４４で算出された結合度ベクトル
【数１２１】
2007043462.tif 000123
【０２４６】
によって該基準値を修正し、修正した基準値に基づき、前処理後の脳機能データρ(l,m,k)に対してボクセル毎に該データ検定を行う。【０２４７】
図２１は、図２０に示した脳機能解析装置４０によって行われるデータ解析の処理手順の一例を示したフローチャートである。【０２４８】
ステップＳ３００において、脳機能データ取得手段１は、MRI装置５０で計測された脳機能情報の原時系列データρ’(l,m,k,i)を取得する。同ステップにおいて、拡散テンソル情報取得手段２は、同じくMRI装置５０で計測された拡散テンソルデータD(l,m,k)を取得する。【０２４９】
次に、ステップＳ３１０において、データ前処理手段３は、ステップＳ３００において取得した脳機能情報の原時系列データρ’(l,m,k,i)に対して、第１の実施形態で詳述したような前処理（ａ）〜（ｄ）を施すことによって、脳機能情報の時系列データρ(l,m,k,i)を生成する。【０２５０】
次に、ステップＳ３２０において、ボクセル間結合度算出手段４４は、ステップＳ３１０において前処理された脳機能情報の時系列データρ(l,m,k,i)をボクセル毎に、例えば、ｚ検定する際の基準値（検定値や有意水準など）を調整するための結合度ベクトル
【数１２２】
2007043462.tif 000124
【０２５１】
を、式（１２）の平均拡散度ベクトル
【数１２３】
2007043462.tif 000125
【０２５２】
の関数として、ステップＳ３１０で取得した拡散テンソルデータD(l,m,k)から算出する。【０２５３】
次に、ステップＳ３３０において、データ検定手段４００は、ステップＳ３１０において前処理された脳機能情報の時系列データρ(l,m,k,i)をボクセル毎にｚ検定を行いｚスコアを算出する。そして、同ステップにおいて、算出したｚスコアが大きなボクセル（例えば、ボクセル(l,m,k)）を基点にして、その周辺（２６方位）のボクセル（例えば、ボクセル(l+1,m,k)）における有意水準の値a(l+1,m,k)を、ステップＳ３２０で算出した結合度ベクトル
【数１２４】
2007043462.tif 000126
【０２５４】
に応じて、例えば、
【数１２５】
2007043462.tif 000127
【０２５５】
のように下げる。ここで、αは重み係数である。【０２５６】
次に、ステップＳ３４０において、画像生成手段７は、ステップＳ３３０において修正された有意水準に基づき、有意であると判断されたボクセル上の脳機能データρ(l,m,k)に基づく値を画像データとして生成する。【０２５７】
そして、ステップＳ３５０において、画像表示手段８は、ステップＳ３４０において生成された画像データを立体表示する。【０２５８】
なお、上記においては、拡散テンソルデータ取得手段２で取得した拡散テンソルデータD(l,m,k)から算出した式（３８）の結合度ベクトル
【数１２６】
2007043462.tif 000128
【０２５９】
を用いて検定の有意水準の値を調整したが、その代わりに、検定値を調整してもよい。【０２６０】
例えば、タスク画像とレスト画像との間でボクセル毎にz検定を行う場合を考えてみる。タスク画像とは、あるタスクを行っているときの画像（例えば、図２の太線で表された１〜３のスライス画像）のことをいい、レスト画像とは、該タスクを行っていないときの画像（例えば、図２の太線で表された５〜７のスライス画像）のことをいう。このとき、ボクセル(l,m,k)におけるzスコア（z検定の検定値）は、通常、
【数１２７】
2007043462.tif 000129
【０２６１】
と表される。ここで、Mt(l,m,k)は、タスク画像の脳機能情報の時系列データρ(l,m,k,i)の平均値、Mc(l,m,k)は、レスト画像の脳機能情報の時系列データρ(l,m,k,i)の平均値、σt(l,m,k)はタスク画像の脳機能情報の時系列データρ(l,m,k,i)の標準偏差、σc(l,m,k)は、レスト画像の脳機能情報の時系列データρ(l,m,k,i)の標準偏差を表す。この前提の下、z=0ならば、タスク画像における脳機能情報の時系列データρ(l,m,k,i)の平均値とレスト画像の脳機能情報の時系列データρ(l,m,k,i)の平均値との間に差がないといえる。Z=1ならば、タスク画像における脳機能情報の時系列データρ(l,m,k,i)の平均値とレスト画像の脳機能情報の時系列データρ(l,m,k,i)の平均値との差が時系列データのばらつきの分だけ大きいことを意味する。このような傾向から、zスコアの値が大きいボクセルがタスクに関連している活性部位であると考えることができる。【０２６２】
そこで、式（４０）を
【数１２８】
2007043462.tif 000130
【０２６３】
のように変更して、ステップＳ３３０において、b(l,m,k)を、ステップＳ３２０で算出した式（３８）の結合度ベクトル
【数１２９】
2007043462.tif 000131
【０２６４】
を用いて、例えば、
【数１３０】
2007043462.tif 000132
【０２６５】
と修正する。ここで、βは重み係数である。【０２６６】
次いで、ステップＳ３４０において、画像生成手段７は、全ボクセルにおいて共通の有意水準に基づき、有意であると判断されたボクセル上の脳機能データρ(l,m,k)に基づく値を画像データとして生成する。そして、ステップＳ３５０において、画像表示手段８は、ステップＳ３４０において生成された画像データを立体表示する。【０２６７】
このように、第４の実施形態においては、MRI装置から取得した拡散テンソルデータに基づいて、同じくMRI装置から取得した脳機能情報の時系列データを検定する際の基準値（検定値又は有意水準）の値をボクセル毎に調整したことによって、脳の活性部位間の結合構造をも考慮した脳の活性部位の特定が可能となる。【０２６８】
＜第４の実施形態の変形例＞
第４の実施形態では、脳機能情報の時系列データρ(l,m,k,i)をボクセル毎に各種データ検定（例えば、Ｔ検定、Ｆ検定、又はＺ検定）する際に、拡散テンソルデータ取得手段２で取得した拡散テンソルデータD(l,m,k)から算出した結合度ベクトル
【数１３１】
2007043462.tif 000133
【０２６９】
を用いて該データ検定の基準値（検定値や有意水準など）を調整したわけだが、本発明の目的を達成するにあたって、脳機能情報の時系列データρ(l,m,k,i)をボクセル毎に各種の解析手法でデータ解析したときに、該解析手法が統計的に有意であるか否かをデータ検定（例えば、Ｔ検定、Ｆ検定、又はＺ検定）するような状況も考えられる。そこで、本変形例として、その際の基準値（検定値や有意水準など）を、結合度ベクトル
【数１３２】
2007043462.tif 000134
【０２７０】
を用いて調整することを考える。なお、重複を避けるため、以降においては第４の実施形態と異なる部分についてのみ説明を加える。【０２７１】
図２２は、第４の実施形態の変形例に係る脳機能解析装置の概略構成を示した概念図である。本変形例に係る脳機能解析装置４０’の構成は、第４の実施形態の構成のデータ検定手段４００の代わりに、データ検定手段４００’とデータ解析手段４６’とを加えた構成を成す。【０２７２】
データ解析手段４６’は、脳機能データ取得手段１で取得され、データ前処理手段３で前処理された脳機能情報の時系列データρ(l,m,k,i)を、従来の分類、回帰、相関等の解析手法を用いてボクセル毎にデータ解析する。【０２７３】
データ検定手段４００’は、データ解析手段４６’で行われた解析手法が統計的に有意であるか否かを各種データ検定（例えば、Ｔ検定、Ｆ検定、又はＺ検定）に基づき検定する際の基準値（検定値や有意水準など）を、結合度ベクトル
【数１３３】
2007043462.tif 000135
【０２７４】
を用いて調整する。【０２７５】
図２３は、図２２に示した脳機能解析装置によって行われるデータ解析の処理手順の一例を示したフローチャ−トである。【０２７６】
ステップＳ３６０において、データ解析手段４６’は、脳機能データ取得手段１で取得され、データ前処理手段３で前処理された脳機能情報の時系列データρ(l,m,k,i)を、従来の分類、回帰、相関等の解析手法を用いてボクセル毎にデータ解析する。【０２７７】
１．分類的解析手法（判別分析、決定木、サポートベクターマシン等）
１−１：判別分析の場合
説明変数を脳機能情報の時系列データρ(l,m,k,i)、外的基準をタスクとレストの時系列データとして、各種の判別関数を用いてデータ分類する。【０２７８】
１−２：決定木の場合
説明変数（属性）をタスクとレストの時系列データ、外的基準として脳機能情報の時系列データρ(l,m,k,i)を離散化して幾つかのクラスを作り、決定木を用いてデータ分類する。【０２７９】
２．回帰的解析手法（各種の関数（線形、非線形）を用いた回帰分析）
２−１：線形回帰分析の場合
説明変数をタスクとレスト、被説明変数を脳機能情報の時系列データρ(l,m,k,i)として、線形回帰式を用いてデータ解析する。【０２８０】
３．相関的解析手法（推定モデル（関数）と実際のデータとの相関を取る手法）
予め設定した、モデル（例えば、刺激を入力とし、脳機能情報の時系列データを出力とするような関数）と実際の脳機能情報の時系列データρ(l,m,k,i)との相関を取ることによって、データ解析する。【０２８１】
ついで、ステップＳ３２０（第４の実施形態と同様）を経て、ステップＳ３３１において、データ検定手段４００’は、ステップＳ３６０において得られた結果が統計的に有意であるか否かを検定する。このとき、検定の基準値（検定値や有意水準など）を結合度ベクトル
【数１３４】
2007043462.tif 000136
【０２８２】
（式（３８））を用いて調整する。例えば、分類的解析手法で得られた結果に対してf検定を行う場合には、ステップＳ３６０において得られた結果に対してボクセル毎にf検定を行いf値を算出する。そして、同ステップにおいて、算出したf値が大きなボクセル（例えば、ボクセル(l,m,k)）を基点にして、その周辺（２６方位）のボクセル（例えば、ボクセル(l+1,m,k)）における有意水準の値a(l+1,m,k)を、ステップＳ３２０で算出した結合度ベクトル
【数１３５】
2007043462.tif 000137
【０２８３】
に応じて、例えば、式（３９）を用いて調整する。回帰的解析手法の場合には、閾値（有意水準に相当）の値を上記のようにして調整すればよいし、相関的解析手法の場合には、相関の閾値を上記のようにして調整すればよい。【０２８４】
次に、ステップＳ３４０において、画像生成手段７は、ステップＳ３３０において修正された有意水準に基づき、有意であると判断されたボクセル上の脳機能データρ(l,m,k)に基づく値を画像データとして生成する。【０２８５】
そして、ステップＳ３５０において、画像表示手段８は、ステップＳ３４０において生成された画像データを立体表示する。【０２８６】
このようにして、本変形例においても、第３の実施形態と同様の効果が得られる。【０２８７】
以上で説明した実施形態はあくまでも例示的なものであり、これに限定されるものではない。【０２８８】
例えば、各実施形態において、隣接ボクセル間の結合度を算出するステップを脳機能データ及び拡散テンソルデータを取得するステップのすぐ後に持ってきてもよい。【０２８９】
また、データ解析の手法についても、ウェブレット解析、ロジスティック回帰分析解析等を適用することも可能である。【０２９０】
また、第１の実施形態に第２又は第３の実施形態を組み込むことも可能である。まず、第１の実施形態に第２の実施形態を組み込む場合には、図１に示した脳機能解析装置１０に、図１０又は図１３に示したデータ平滑化手段２００又は２００’と同じ機能を有するデータ平滑化手段を更に備えるような構成とすればよい。そして、第１の実施形態のステップＳ５０において、データ解析手段６がデータ評価量Qを最小にするような回帰係数を求めた後に、該回帰係数の値に対して上記データ平滑化手段を適用することにより該回帰係数の値を異方的に平滑化すればよい。これにより、第１の実施形態で得られる画像と比較した場合、白質（神経繊維）の活性度がより強調された画像を生成することが可能となる。【０２９１】
また、第１の実施形態に第３の実施形態を組み込む場合には、図１に示した脳機能解析装置１０に、図１５又は図１８に示したクラスタリング手段３００又は３００‘と同じ機能を有するクラスタリング手段を更に備えるような構成とすればよい。そして、第１の実施形態のステップＳ５０において、データ解析手段６がデータ評価量Qを最小にするような回帰係数を求めた後に、該回帰係数の値に対して上記クラスタリング手段を適用することによりボクセル同士をクラスタリングさせればよい。これにより、第１の実施形態で得られる画像と比較した場合、白質（神経線維）の活性度がより強調された画像を生成することが可能となる。【０２９２】
さらに、上記した各実施形態で得られた結果に対して、さらに、拡散テンソルデータD(l,m,k)を用いて、白質ボクセルの活性化部位を伸ばして行くことも可能である。その手法は基本的にトラクトグラフィー（tractography）と同様の手法である。【０２９３】
さらに、上記の手法（回帰分析、相関、分類、独立成分分析、検定等）を白質ボクセルだけに適用することで、白質の活性化を捕らえることも可能である。【０２９４】
さらに、各実施形態における実験としてはブロックデザインを想定したが、本発明は、事象関連ｆMRIにも適用可能である。【０２９５】
また、本発明は接続解析（Connectivity Analysis）にも適用可能である。接続解析の手法にはいくつかあるが、たとえば、SEM(Structural Equation Modeling)(あるいは共分散構造分析)の場合には、その計算過程で相関や検定が行われるが、その際に本発明を適用できるし、平滑化やクラスタリングも適用できる。【０２９６】
上記では、神経細胞と神経線維の結合に関する情報を拡散テンソルデータから獲得したが、すでに解剖学的に既知となっている結合に関する情報（知見）があれば、それを用いることも可能である。【０２９７】
本発明の趣旨を逸脱しない限りにおいて、その他にも様々な実施形態が可能であり、本発明はあくまでも特許請求の範囲によってのみ規定されるものである。【０２９８】
最後に、本発明の脳機能解析方法を用いた脳機能解析の効果を示す。図２４は、被験者に聴覚刺激を用い単純復唱を行わせた場合の各種脳機能解析の結果を示したものである。実験は、ヘッドホンから流れてくる音声を聞こえた通りに復唱することをタスクとし、レストは出来るだけ何も考えない状態とした。タスクのボリューム数（計測回数）は48，レストのボリューム数（計測回数）は60，解析に使用しないボリューム数（計測回数）は12とした（図２を参照）。また、TR（Time of Repetition：図２のボリューム幅に対応）は5000msとした。被験者の数は８人であるが、SPMと本発明の脳機能解析方法の両方で良好な結果を得た被験者の解析結果を以下に示す。【０２９９】
図２４Ａ，２４Ｂは共にSPMを用いたT検定の結果を表すものであり、図２４ＡはT検定の閾値に修正を施した場合の解析結果を表し、図１２ＢはT検定の閾値に修正を施さなかった場合の解析結果を表している。図２４Ｃは、SPMとトラクトグラフィーとを組み合わせた解析手法による解析結果を表し、図２４Ｄは本発明の第１の実施形態に係る脳機能解析方法を用いた解析結果を表す。なお、SPMの前処理として、位置補正、正規化、平滑化を行った。平滑化の半値幅は9mmとした。また、トラクトグラフィーには拡散テンソル解析用のソフトウェアとしてdTVを用いた。【０３００】
前提として、単純復唱の時には、ウェルニッケ野（感覚性言語野）とブローカー野（運動性言語野）が活性化し、両野は弓状束という神経繊維の束で結合されていることが解剖学的に知られている。この弓状束は発話の際に用いられている。【０３０１】
図２４Ａに示すように、SPMを用いた場合でも、T検定の閾値を修正した場合にはウェルニッケ野、ブローカー野の活性化は検出できない。一方、図２４Ｂに示すように、SPMを用い、T検定の閾値を修正しなかった場合には、両野の活性化は検出できたが、両野をつなぐ弓状束は検出できなかった。これらの結果に対し、図２４Ｃに示すように、SPMとトラクトグラフィーとを組み合わせた解析手法を用いると、ウェルニッケ野及びブローカー野の活性化の検出に加え、弓状束を検出できた。【０３０２】
一方、本発明の第１の実施形態に係る脳機能解析方法を用いた場合には、図２４Ｄに示すように、両野の活性化や両野をつなぐ弓状束の検出に加えて、弓状束の活性化（弓状束を構成する神経線維におけるパルス伝達）まで検出できる。【０３０３】
また、SPMとトラクトグラフィーとを組み合わせた解析手法においては、神経繊維（この場合、弓状束）の検出にあたり、解析者が神経繊維の追跡開始点と終了点とを予め指定しなければならない。したがって、神経線維の追跡開始点と終了点とが予め決定できないような場合（例えば、灰白質の活性領域が予めわかっていないようなタスク）にはこの手法は有効でないことになる。それに対し、本発明の第１の実施形態に係る脳機能解析方法によれば、神経線維の追跡開始点と終了点とが予め指定できない場合に対しても、図２４Ｄに示したように、神経線維の活性化の検出が可能である。
【産業上の利用可能性】
【０３０４】
このように、本発明によれば、脳の活性部位を特定できるのはもちろんのこと、解析者が予め神経線維の開始点及び終了点を指定しなくとも、神経線維（活性部位間の結合構造）を特定することが可能となり、さらには神経繊維の活性化（神経線維におけるパルス伝達）をも検出可能となるので、認知症、失語症、統合失調症等の高次脳機能障害の診断及び治療に対して極めて有効な脳機能解析方法、及び脳機能解析プログラムを提供することができる。
【図１】
2007043462.tif 000138
【図２】
2007043462.tif 000139
【図３】
2007043462.tif 000140
【図４】
2007043462.tif 000141
【図５】
2007043462.tif 000142
【図６】
2007043462.tif 000143
【図７】
2007043462.tif 000144
【図８】
2007043462.tif 000145
【図９】
2007043462.tif 000146
【図１０】
2007043462.tif 000147
【図１１】
2007043462.tif 000148
【図１２】
2007043462.tif 000149
【図１３】
2007043462.tif 000150
【図１４】
2007043462.tif 000151
【図１５】
2007043462.tif 000152
【図１６】
2007043462.tif 000153
【図１７】
2007043462.tif 000154
【図１８】
2007043462.tif 000155
【図１９】
2007043462.tif 000156
【図２０】
2007043462.tif 000157
【図２１】
2007043462.tif 000158
【図２２】
2007043462.tif 000159
【図２３】
2007043462.tif 000160
【図２４】
2007043462.tif 000161
【手続補正書】
【提出日】平成20年11月13日(2008.11.13)
【手続補正１】
【補正対象書類名】特許請求の範囲
【補正対象項目名】全文
【補正方法】変更
【補正の内容】
【特許請求の範囲】
【請求項１】
脳の活性部位を特定可能とする脳機能データと前記脳内のプロトンの拡散度を特定可能とする拡散テンソルデータとをボクセル単位で取得する第１のステップと、
前記拡散テンソルデータに基づき隣接ボクセル間の結合度の評価量を構成する第２のステップと、
前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データの解析を行う第３のステップと、
を有することを特徴とする脳機能データ解析方法。
【請求項２】
請求項１に記載の脳機能解析方法において、前記第３のステップは、前記脳機能データ及びタスクの一方を説明変数とし、前記脳機能データ及び前記タスクの他方を被説明変数とする回帰分析であって、前記脳機能データのボクセル毎の評価量に前記隣接ボクセル間の結合度の評価量を組み込んだ評価量の最小値及び最大値の一方を算出することを特徴とする脳機能データ解析方法。
【請求項３】
請求項１に記載の脳機能解析方法において、前記第３のステップは、前記脳機能データを検定する手法であって、前記検定の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能データ解析方法。
【請求項４】
請求項１に記載の脳機能解析方法において、前記第３のステップは、前記脳機能データを分類する手法であって、前記分類の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能データ解析方法。
【請求項５】
請求項１に記載の脳機能解析方法において、前記第３のステップは、前記脳機能データと所定のモデルとの相関をとる手法であって、前記相関の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能データ解析方法。
【請求項６】
請求項１に記載の脳機能解析方法において、前記第３のステップは、前記脳機能データから主要な成分を抽出する手法であって、前記抽出の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能データ解析方法。
【請求項７】
請求項１に記載の脳機能解析方法において、前記第３のステップにおいて、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データを平滑化することを特徴とする脳機能データ解析方法。
【請求項８】
請求項１に記載の脳機能解析方法において、前記第３のステップにおいて、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データをクラスタリングすることを特徴とする脳機能データ解析方法。
【請求項９】
請求項１乃至８のいずれかに記載の脳機能解析方法において、前記第１のステップにおいて取得した拡散テンソルデータに基づき、同ステップにおいて取得した脳機能データから特定される活性部位ボクセルの脳機能データの値を他のボクセルの脳機能データの値に伝播させる前処理を施す第４のステップをさらに備えたことを特徴とする脳機能データ解析方法。
【請求項１０】
コンピュータを、
脳の活性部位を特定可能とする脳機能データをボクセル単位で取得する脳機能データ取得手段と、
前記脳内のプロトンの拡散度を特定可能とする拡散テンソルデータをボクセル単位で取得する拡散テンソルデータ取得手段と、
前記拡散テンソルデータに基づき隣接ボクセル間の結合度の評価量を構成するデータ評価量構成手段と、
前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データの解析を行うデータ解析手段と、
して機能させることを特徴とする脳機能解析プログラム。
【請求項１１】
請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記脳機能データ及びタスクの一方を説明変数とし、前記脳機能データ及び前記タスクの他方を被説明変数とする回帰分析する手段であって、前記脳機能データのボクセル毎の評価量に前記隣接ボクセル間の結合度の評価量を組み込んだ評価量の最小値及び最大値の一方を算出することを特徴とする脳機能解析プログラム。
【請求項１２】
請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記脳機能データを検定する手段であって、前記検定の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析プログラム。
【請求項１３】
請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記脳機能データを分類する手段であって、前記分類の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析プログラム。
【請求項１４】
請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記脳機能データと所定のモデルとの相関をとる手段であって、前記相関の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析プログラム。
【請求項１５】
請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記脳機能データから主要な成分を抽出する手段であって、前記抽出の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析プログラム。
【請求項１６】
請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析は、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データを平滑化することを特徴とする脳機能解析プログラム。
【請求項１７】
請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データをクラスタリングすることを特徴とする脳機能解析プログラム。
【請求項１８】
請求項１０乃至１７のいずれかに記載の脳機能解析プログラムにおいて、前記コンピュータを、前記拡散テンソルデータ取得手段によって取得した拡散テンソルデータに基づき、前記脳機能データ取得手段によって取得した脳機能データから特定される活性部位ボクセルの脳機能データの値を他のボクセルの脳機能データの値に伝播させる前処理を施すデータ前処理手段として機能させることを特徴とする脳機能解析プログラム。
【請求項１９】
脳の活性部位を特定可能とする脳機能データをボクセル単位で取得する脳機能データ取得手段と、
前記脳内のプロトンの拡散度を特定可能とする拡散テンソルデータをボクセル単位で取得する拡散テンソルデータ取得手段と、
前記拡散テンソルデータに基づき隣接ボクセル間の結合度の評価量を構成するデータ評価量構成手段と、
前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データの解析を行うデータ解析手段と、
を備えたことを特徴とする脳機能解析装置。
【請求項２０】
請求項１９に記載の脳機能解析装置において、前記データ解析手段は、前記脳機能データ及びタスクの一方を説明変数とし、前記脳機能データ及び前記タスクの他方を被説明変数とする回帰分析する手段であって、前記脳機能データのボクセル毎の評価量に前記隣接ボクセル間の結合度の評価量を組み込んだ評価量の最小値及び最大値の一方を算出することを特徴とする脳機能解析装置。
【請求項２１】
請求項１９に記載の脳機能解析装置において、前記データ解析手段は、前記脳機能データを検定する手段であって、前記検定の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析装置。
【請求項２２】
請求項１９に記載の脳機能解析装置において、前記データ解析手段は、前記脳機能データを分類する手段であって、前記分類の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析装置。
【請求項２３】
請求項１９に記載の脳機能解析装置において、前記データ解析手段は、前記脳機能データと所定のモデルとの相関をとる手段であって、前記相関の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析装置。
【請求項２４】
請求項１９に記載の脳機能解析装置において、前記データ解析手段は、前記脳機能データから主要な成分を抽出する手段であって、前記抽出の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析装置。
【請求項２５】
請求項１９に記載の脳機能解析装置において、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データを平滑化するデータ平滑化手段をさらに備えたことを特徴とする脳機能解析装置。
【請求項２６】
請求項１９に記載の脳機能解析装置において、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データをクラスタリングするデータクラスタリング手段をさらに備えたことを特徴とする脳機能解析装置。
【請求項２７】
請求項１９乃至２６のいずれかに記載の脳機能解析装置において、前記拡散テンソルデータ取得手段によって取得した拡散テンソルデータに基づき、前記脳機能データ取得手段によって取得した脳機能データから特定される活性部位ボクセルの脳機能データの値を他のボクセルの脳機能データの値に伝播させる前処理を施すデータ前処理手段をさらに備えたことを特徴とする脳機能解析装置。
【手続補正３】
【補正対象書類名】明細書
【補正対象項目名】０００１
【補正方法】変更
【補正の内容】
【０００１】
本発明は、fMRI（functional Magnetic Resonance Imaging：機能的磁気共鳴撮影）やPET（Positron Emission Tomography：陽電子断層撮影）等の計測手法により得られる各種のデータを用いて脳機能を解析するための脳機能データ解析方法、脳機能解析装置及び脳機能解析プログラムに関する。
【手続補正４】
【補正対象書類名】明細書
【補正対象項目名】００１５
【補正方法】変更
【補正の内容】
【００１５】
本発明の一つの側面は、脳の活性部位を特定可能とする脳機能データと前記脳内のプロトンの拡散度を特定可能とする拡散テンソルデータとをボクセル単位で取得する第１のステップと、前記拡散テンソルデータに基づき隣接ボクセル間の結合度の評価量を構成する第２のステップと、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データの解析を行う第３のステップとを有することを特徴とする脳機能データ解析方法を提供することである。
本発明の他の側面は、脳の活性部位を特定可能とする脳機能データをボクセル単位で取得する脳機能データ取得手段と、前記脳内のプロトンの拡散度を特定可能とする拡散テンソルデータをボクセル単位で取得する拡散テンソルデータ取得手段と、前記拡散テンソルデータに基づき隣接ボクセル間の結合度の評価量を構成するデータ評価量構成手段と、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データの解析を行うデータ解析手段とを備えたことを特徴とする脳機能解析装置を提供することである。
【手続補正５】
【補正対象書類名】明細書
【補正対象項目名】０１１７
【補正方法】変更
【補正の内容】
【０１１７】
のように求めることができる。ここで、Xはρ(n,i)を成分とするL3×I行列
【数６２】
2007043462.tif 000162
【手続補正６】
【補正対象書類名】明細書
【補正対象項目名】０１４２
【補正方法】変更
【補正の内容】
【０１４２】
＜第１の実施形態の変形例４＞
第１の実施形態では、ステップＳ５０において、データ解析手段６は、式（２４）の線形回帰式を用いることによりノンパラメトリック回帰分析をしたわけだが、データ解析の手法として、判別分析（非特許文献８：“多変量解析のはなし”，有馬哲，石村貞夫， 東京図書，1987）、数量化II類（非特許文献９：“多変量解析のはなし”，有馬哲，石村貞夫， 東京図書，1987）、決定木（非特許文献１０：“AIによるデータ解析”，J.R.キンラン，トッパン，1995）、サポートベクターマシン（非特許文献１１：“サポートベクターマシン入門”，Nello Cristianini and John Shawe‐Taylor，共立出版，2005）等の分類的解析手法を用いることもできる。
【手続補正７】
【補正対象書類名】明細書
【補正対象項目名】０２９９
【補正方法】変更
【補正の内容】
【０２９９】
図２４Ａ，２４Ｂは共にSPMを用いたT検定の結果を表すものであり、図２４ＡはT検定の閾値に修正を施した場合の解析結果を表し、図２４ＢはT検定の閾値に修正を施さなかった場合の解析結果を表している。図２４Ｃは、SPMとトラクトグラフィーとを組み合わせた解析手法による解析結果を表し、図２４Ｄは本発明の第１の実施形態に係る脳機能解析方法を用いた解析結果を表す。なお、SPMの前処理として、位置補正、正規化、平滑化を行った。平滑化の半値幅は9mmとした。また、トラクトグラフィーには拡散テンソル解析用のソフトウェアとしてdTVを用いた。
【国際調査報告】
2007043462.tif 000163
2007043462.tif 000164
2007043462.tif 000165
2007043462.tif 000166
2007043462.tif 000167
2007043462.tif 000168
2007043462.tif 000169