JP4308871B2 - 脳機能データ解析方法、脳機能解析装置及び脳機能解析プログラム - Google Patents

脳機能データ解析方法、脳機能解析装置及び脳機能解析プログラム Download PDF

Info

Publication number
JP4308871B2
JP4308871B2 JP2007539917A JP2007539917A JP4308871B2 JP 4308871 B2 JP4308871 B2 JP 4308871B2 JP 2007539917 A JP2007539917 A JP 2007539917A JP 2007539917 A JP2007539917 A JP 2007539917A JP 4308871 B2 JP4308871 B2 JP 4308871B2
Authority
JP
Japan
Prior art keywords
brain function
data
analysis
degree
brain
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
JP2007539917A
Other languages
English (en)
Other versions
JPWO2007043462A1 (ja
Inventor
洋 月本
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tokyo Denki University
Original Assignee
Tokyo Denki University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Tokyo Denki University filed Critical Tokyo Denki University
Publication of JPWO2007043462A1 publication Critical patent/JPWO2007043462A1/ja
Application granted granted Critical
Publication of JP4308871B2 publication Critical patent/JP4308871B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00—Arrangements or instruments for measuring magnetic variables
    • G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
    • G01R33/44—Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
    • G01R33/48—NMR imaging systems
    • G01R33/54—Signal processing systems, e.g. using pulse sequences ; Generation or control of pulse sequences; Operator console
    • G01R33/56—Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution
    • G01R33/563—Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution of moving material, e.g. flow contrast angiography
    • G01R33/56341—Diffusion imaging
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/05—Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves 
    • A61B5/055—Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves  involving electronic [EMR] or nuclear [NMR] magnetic resonance, e.g. magnetic resonance imaging
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/40—Detecting, measuring or recording for evaluating the nervous system
    • A61B5/4058—Detecting, measuring or recording for evaluating the nervous system for evaluating the central nervous system
    • A61B5/4064—Evaluating the brain
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis
    • A61B5/7264—Classification of physiological signals or data, e.g. using neural networks, statistical classifiers, expert systems or fuzzy systems
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00—Arrangements or instruments for measuring magnetic variables
    • G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
    • G01R33/44—Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
    • G01R33/48—NMR imaging systems
    • G01R33/4808—Multimodal MR, e.g. MR combined with positron emission tomography [PET], MR combined with ultrasound or MR combined with computed tomography [CT]
    • G01R33/481—MR combined with positron emission tomography [PET] or single photon emission computed tomography [SPECT]
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00—Arrangements or instruments for measuring magnetic variables
    • G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
    • G01R33/44—Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
    • G01R33/48—NMR imaging systems
    • G01R33/54—Signal processing systems, e.g. using pulse sequences ; Generation or control of pulse sequences; Operator console
    • G01R33/56—Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution
    • G01R33/5608—Data processing and visualization specially adapted for MR, e.g. for feature analysis and pattern recognition on the basis of measured MR data, segmentation of measured MR data, edge contour detection on the basis of measured MR data, for enhancing measured MR data in terms of signal-to-noise ratio by means of noise filtering or apodization, for enhancing measured MR data in terms of resolution by means for deblurring, windowing, zero filling, or generation of gray-scaled images, colour-coded images or images displaying vectors instead of pixels

Description

本発明は、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:傾斜磁場）：
を印加することによって拡散の異方性を計測する手法であり、例えば、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)は、
と書くことができ、拡散強調画像を生成する際には式（１）における拡散テンソル
がデータ解析の対象となる。ここで、ρ’(l,m,k,i)はMPGを印加しない場合のBOLD信号(通常の脳機能解析におけるデータ解析の対象)の強度を表し、bはMPGの強さを表すパラメータである。なお、ρ’(l,m,k,i)は、
と表される。ここで、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)から隣接ボクセル間の結合度を表す結合度ベクトル
を算出するボクセル間結合度算出手段４と、データ前処理手段３で前処理された脳機能情報の時系列データρ(l,m,k,i)とボクセル間結合度算出手段４で算出された結合度ベクトル
とから所定の評価量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)は式（２）に示したように行列で表現され、例えば、
のような値をとる。この例で示すように、拡散テンソルは通常、対称テンソル（非対角成分が対称）なので、実質的には６個の成分を有することになる。また、図中の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)の固有方程式
を満たす固有ベクトル
の中で最大固有値λMに対応する固有ベクトル
を算出する（α=1,2,3）。説明を簡単にするために、以降においては、｜λ1｜≧｜λ2｜≧｜λ3｜とし、最大固有値であるλ1に対応する固有ベクトル
を
と書くことにする。また、該固有ベクトル
のノルムは１に規格化されているものとする。
次いで、同ステップにおいて、データ前処理手段３は、まず、あるボクセル(l,m,k)に対して、該ボクセル(l,m,k)と隣接する２６個のボクセル（３次元）の中から該ボクセル(l,m,k)上の固有ベクトル
との内積の絶対値が最大となる方向ベクトルに対応する隣接ボクセルを求める。この条件を満たす隣接ボクセルを、例えば、(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)上の固有ベクトル
との内積が最大となるボクセルであるから、ボクセル(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)の値を、
と書き換える。
（状況２）予め取得された所定の脳の解剖学的データ（MRI装置５０により取得されたT1強調画像）から、ボクセル(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)が共に灰白質とは見なされない（つまり、白質（神経の神経繊維：図中のW）と見なされる）場合（図９（Ｃ）を参照）：
この場合は、ボクセル(l,m,k)上の脳機能情報の時系列データρ’(l,m,k,i)の値を、
と書き換える。
このような処理を全ボクセルに対して行う。その結果、図９（Ａ），（Ｂ）の下図に示したように、灰白質と見なされるボクセル（灰白質ボクセル）上の脳機能情報の時系列データの値を白質と見なされる隣接ボクセル（白質ボクセル）上の脳機能情報の時系列データの値に反映（伝播）させることができる。また、図９（Ｃ）の下図に示したように、白質と見なされるボクセル同士が隣接している場合には、これらのボクセル上の脳機能データの値が平滑化されることになる。このような前処理を施すことによって、以降のステップの効果を際立たせることができる。
なお、このような前処理（ｄ）の仕方は、上記した例に限らず、例えば、灰白質ボクセルに隣接し、且つ該灰白質ボクセルの脳機能情報の時系列データの値を受け取った白質ボクセルが、さらに、隣接する白質ボクセルに対して上記と同様の手順によってその値を順次伝播させるようにしてもよい。また、上記のように灰白質ボクセルの脳機能情報の時系列データの値を白質ボクセルに渡す際に、拡散テンソルデータD(l,m,k)に応じて重み付けをした値を渡すような構成にしてもよい。さらに、上記の手順で書き換えられた白質ボクセルの値と元々の白質ボクセルの値のうちの大きい方を白質ボクセルの値としてもよい。
次に、ステップＳ３０において、ボクセル間結合度算出手段４は、上記の拡散テンソルデータ
に基づき、隣接ボクセル間の結合度ベクトル
を、例えば、
から算出する。ここで、
の各成分は、例えば、
と定義されるのものであり、以降において、これを隣接ボクセル間の平均拡散度ベクトルと称する。αは
を満たすパラメータであり、
は
を満たす定数ベクトルである。
式（１３）中の
は、プロトンの拡散度を表すベクトル（以降、これを拡散度ベクトルと称する）であり、第１成分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方向とは図７に示した方向と同一の方向を想定している。なお、式（９）においては、拡散テンソルが対称テンソルであることを前提としている。
このような意味合いを有する拡散度ベクトル
を用いて、式（１３）を介して式（１１）のように構成された隣接ボクセル間の結合度ベクトル
の各成分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)）を介して間接的に考慮されている。もちろん、斜め方向に隣接するボクセル間の結合度を直接考慮するような構成にしてもよい。
さらに、拡散度ベクトル
としては、例えば、
（Ｉ）楕円体モデルにおける
を採用する。ここで、｜D｜は、式（９）で表される行列の行列式を表す。また、例えば、
（II）式（９）で表される行列の対角成分Dll(l,m,k)，Dmm(l,m,k)，Dkk(l,m,k)を用いて、
とすることもできる。また、例えば、
（III）FA（Fractional Anisotropy）：
やRA（Relational Anisotropy）：
に代表される拡散の異方性を表す指標を用いて拡散度ベクトル
を決めることも可能である。ここで、DAVは、
である。その他にも、
（IV）予め取得された所定の脳の解剖学的データ（MRI装置５０により取得されたT1強調画像）から脳の白質ボクセルと灰白質ボクセルとを判断し、白質ボクセルでは上記（Ｉ）の拡散度ベクトルを用い、灰白質ボクセルでは、その他の拡散度ベクトル（例えば、上記（II）又は（III））を用いることも可能である。
このように、本発明における拡散度ベクトル
としては、MRI装置５０から取得可能な拡散テンソルデータD(l,m,k)から神経線維の走行方向が特定可能な量であればいかなるものを用いてもよい。
また、パラメータαに関する不等式（１４）の代わりに、例えば、灰白質ボクセルで大きく、白質ボクセルで小さい値といったようにパラメータαの値を灰白質ボクセルと白質ボクセルとで切替えるように設定をしてもよい。
また、式（１１）の結合度ベクトル
の代わりに、
を用いてもよい。
さらに、神経繊維の走行方向を強調するために、式（１１）又は式（２２）の結合度ベクトル
の各成分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以上の値を有する成分のみを例えば上記の方法で強調してもよい。
さらに、灰白質ボクセルは隣接する６個のボクセル（６方位）とのみ結合させ、白質ボクセルは固有方程式（５）の最大固有値に対応する固有ベクトルとの内積が最大のボクセル及びその点対称のボクセル（２方位）とのみ結合するように結合度ベクトル
を構成してもよい。
次に、ステップＳ４０において、データ評価量構成手段５は、ステップＳ２０において前処理された脳機能情報の時系列データρ(l,m,k,i)とステップＳ４０において算出された結合度ベクトル
とから、図５（Ａ）に例示したようなデータをノンパラメトリック回帰分析するために、例えば、隣接ボクセルを６個とした場合には、データ評価量Q：
を構成する。ここで、iは時間、l,m,kはそれぞれボクセルのX方向、Y方向及びZ方向の位置を表す変数であることは上述した通りである。ただし、実際には、図７に示す立体構成されたボクセルの中には脳（神経）を表さないボクセルもあるので、その領域に関するl，m，kは除外するものとする。また、φ(i)は、図５（Ａ）におけるCLASSを表す変数の時系列データであり、この場合、1（=ON）と0（=OFF）のニ値をとる。本実施形態に係るノンパラメトリック回帰分析では、図５（Ａ）に示したような脳機能情報の時系列データρ(l,m,k,i)を説明変数、CLASSを表す変数の時系列データφ(i)を被説明変数とし、これらの変数の間に、
を回帰係数とした線形回帰式：
を想定する。
及び
はそれぞれ推定解を表し、κは隣接ボクセル間の連続性を表す重みパラメータである。このように、式（２３）は隣接ボクセル間の連続性に基づいた評価量を構成したが、その代わりに、隣接ボクセル間の滑らかさに基づいた評価量を構成してもよい。
本実施形態におけるデータ解析の手法としてノンパラメトリック回帰分析を用いる根拠は、被説明変数φ(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’：
を最小にするような
を求める。ここで、右辺第１項は、各時刻iにおける被説明変数φ(i)の推定解
からの誤差を表す項であり、第２項が、被説明変数φ(i)の時間方向に対する連続性を評価するための項である。しかしながら、図５（Ａ）に示す例からわかるように被被説明変数φ(i)の値に連続性を仮定することはできない。そこで、本実施形態では、このような連続性の制約を被説明変数φ(i)に対してではなく、説明変数ρ(l,m,k,i)に対して課すように拡張したノンパラメトリック回帰分析（非特許文献７：“脳機能画像のノンパラメトリック回帰分析”,月本 洋 他,電子情報通信学会誌D-II，Vol. J84, No.12，pp2623-2633，2001年12月）に基づいて、式（２３）で表されるデータ評価量Qを構成する。
ここで、式（２３）の意味について説明する。右辺第１項は、式（２５）と同じく、各時刻iにおける被説明変数φ(i)の推定解
からの誤差を表す項である。第２項〜第４項は、隣接するボクセルが被説明変数φ(i)に与える影響の連続性を表す項であり、この連続性は隣接する説明変数ρ(l,m,k,i)が被説明変数に与える影響、すなわち、式（２４）の回帰係数
の連続性、で評価できる。
さらに、本実施形態においては、第２項〜第４項において、ボクセル間結合度算出手段４で算出された結合度ベクトル
の各成分Cl(l,m,k)，Cm(l,m,k)，Ck(l,m,k)をそれぞれX方向、Y方向、Z方向に対する隣接ボクセル間の連続性に対する重み付け因子として導入している。式（１１）に示したように、これらの成分はそれぞれ、平均拡散度ベクトル
の各成分D’l(l,m,k)，D’m(l,m,k)，D’k(l,m,k)の関数となっているので、結果として、式（２３）のデータ評価量Qは神経線維（神経細胞）の走行方向（プロトンの拡散の異方性）を反映していることになる。
次に、ステップＳ５０において、データ解析手段６は、ステップＳ４０において構成されたデータ評価量Qの極値を求めるための計算を行う。本実施形態においては、式（２３）のデータ評価量Qを最小にするような回帰係数
を算出する。ここでは、簡単のため、式（２３）を
のように書き換えたデータ評価量Q”について説明する。ここで、nは立体構成された全ボクセルに対して所定の番号付けを施すことによって１次元配列し直したボクセルの位置を表す正の離散変数であり、３変数表示(l,m,k)と１対１対応している。さらに、簡単のため、ここではL=M=Kとしたので、全ボクセル数はL3となる。また、C(n)は、このように１次元に配列し直した際に、ボクセルnと該ボクセルnと隣接するボクセルn’との間の結合度を表す係数であり、この係数は、式（１１）の結合ベクトル
から一意に定まる。例えば、図３（Ｂ）に示した２次元スライス画像Skに話を限定した場合、n=130のボクセルと隣接するボクセルはn=127,129,131,257の４つとなる。ステップＳ３０で算出された結合度ベクトル
により、ボクセルn=127とボクセルn=130との結合度、ボクセルn=129とボクセルn=130との結合度、ボクセルn=131とボクセルn=130との結合度、ボクセルn=257とボクセルn=130との結合度はそれぞれ、例えば、2，2，5，5のように決まる。つまり、この場合、
となる。
このとき、回帰係数
は、最小２乗法を用いて、
のように求めることができる。ここで、Xはρ(n,i)を成分とするL3×I行列
はφ(i)を成分とするI次元ベクトル、
は
を成分とするL3次元ベクトル、κOは交差検証法（Cross Validation）で求めたκの最適値、すなわち、
を最小にするようなκである。ここで、
はi番目の事例を除いて回帰係数を計算し、その回帰係数を用いてi番目の事例の測定値を予測した値である。
次に、ステップＳ６０において、画像生成手段７は、上記のようにして算出された回帰係数
を、例えば、回帰係数の値が0（すなわち、無相関）のボクセルをグレイスケールの基準として、回帰係数の値が正（すなわち、正の相関）のボクセルを白方向の色調とし、回帰係数の値が負（すなわち、負の相関）のボクセルを黒方向の色調とするような色調を施した画像を生成する。その他にも、グレイスケール上で取り得るボクセルの値の頻度のヒストグラムにおいて、上位５％の値以上を表示させるようにしてもよいし、所定の規則に従ってフルカラー画像を生成してもよい。
そして、ステップＳ７０において、画像表示手段８は、上記のようにして生成された画像を立体表示する。
このように、第１の実施形態においては、MRI装置から取得した脳機能情報の時系列データに基づく誤差の評価項に加えて、同じくMRI装置から取得した拡散テンソルデータに基づく空間方向の連続性の評価項を加えた評価量を構成した。そして、この評価項をノンパラメトリック線形回帰分析し、その結果得られた回帰係数を画像表示の対象としたことによって、脳の活性部位間の結合構造をも考慮した脳の活性部位の特定が可能となる。
なお、本実施形態においては、ステップＳ２０において通常の前処理（ａ）〜（ｃ）以外に前処理（ｄ）のプロセスを設けたが、前処理（ｄ）をはずしてもよい。このことは、以降の実施形態においても同様である。
なお、第１の実施形態では、データ評価量として式（２３）を構成したわけだが、本発明の目的を達成するにあたり、データ解析の手法に応じて、その他にも各種の評価量を構成することができる。そこで、次に、上記した第１の実施形態の幾つかの変形例について説明するが、重複を避けるため、以降においては第１の実施形態と異なる部分についてのみ説明を加える。
＜第１の実施形態の変形例１＞
第１の実施形態では、ステップＳ５０において、データ解析手段６は、式（１９）の線形回帰式を用いることによりノンパラメトリック回帰分析をしたわけだが、より厳密には、回帰式を２次関数、３次関数等の一般のｎ次関数で行うことも考えられる。このような非線形回帰分析を行う際には、ニューラルネットワーク（Neural Network）モデルを用いることができる。典型的なニューラルネットワークモデルとして、入力層、中間層、出力層から成る階層型ニューラルネットワークモデルがある。
本変形例では、ステップＳ４０において、データ評価量構成手段５は、ニューラルネットワークモデルにおける通常の学習の評価量である２乗誤差の項に対して、ボクセル間結合度算出手段４で算出された結合度ベクトル
に基づく隣接ボクセル間の連続性を評価する評価量を加えた評価量を構成する。そして、ステップＳ５０において、データ解析手段６は、このようにして構成された評価量の最小値を求めるために、誤差逆伝播法（Error Back Propagation Method）等を用いて解析した後に、ステップＳ６０において、画像生成手段７は、入力（脳機能情報の時系列データ）と出力の間の感度（入力を変化させたときに出力に与える影響）を表示画像として第１の実施形態と同様にして生成する。そして、ステップＳ７０において、画像表示手段８は、上記のようにして生成された画像を立体表示する。
これによって、第１の実施形態と同様の効果が得られる。
＜第１の実施形態の変形例２＞
第１の実施形態では、ステップＳ４０において、データ評価量構成手段５は、脳機能の時系列データρ(l,m,k,i)を説明変数、φ(i)を被説明変数として式（２３）の評価量Qを構成したわけだが、本変形例では、その逆、つまり、φ(i)を説明変数、ρ(l,m,k,i)を被説明変数として、ボクセル(l,m,k)毎に、例えば、
といった評価量Q”’を構成する。ここで、
である。
そして、ステップＳ５０において、データ解析手段６は、このようにして構成された評価量Q”’の最小値を求めるために、最適なκを交差検証法を用いて決定し、GA（Genetic Algorithm:遺伝的アルゴリズム）等の確率的探索アルゴリズムによって回帰係数
を求める。
次いで、ステップＳ６０において、画像生成手段７は、上記の解析結果を各種検定した際に有意であると判断されたボクセルの脳機能データρ(l,m,k)に基づく値を表示画像として生成する。そして、ステップＳ７０において、画像表示手段８は、上記のようにして生成された画像を立体表示する。
これによって、第１の実施形態と同様の効果が得られる。
＜第１の実施形態の変形例３＞
第１の実施形態の変形例１では、ステップＳ４０において、データ評価量構成手段５は、脳機能の時系列データρ(l,m,k,i)を説明変数、φ(i)を被説明変数として式（２３）の評価量を構成したわけだが、本変形例では、その逆、つまり、φ(i)を説明変数、ρ(l,m,k,i)を被説明変数として、階層型ニューラルネットワークモデルを用いて、非線形回帰分析を行うこともできる。ただし、この場合には、ステップＳ５０において極値を求める際に誤差逆伝播法が使えないので、その代わりに、データ解析手段６は、最適なκを交差検証法を用いて決定し、GA等の確率的探索アルゴリズムによってデータ解析を行う。そして、ステップＳ６０において、画像生成手段７は、上記の解析結果を各種検定した際に有意であると判断されたボクセルの脳機能データρ(l,m,k)に基づく値を表示画像として生成する。そして、ステップＳ７０において、画像表示手段８は、上記のようにして生成された画像を立体表示する。
これによって、第１の実施形態と同様の効果が得られる。
＜第１の実施形態の変形例４＞
第１の実施形態では、ステップＳ５０において、データ解析手段６は、式（２４）の線形回帰式を用いることによりノンパラメトリック回帰分析をしたわけだが、データ解析の手法として、判別分析（非特許文献８：“多変量解析のはなし”，有馬哲，石村貞夫， 東京図書，1987）、数量化II類（非特許文献９：“多変量解析のはなし”，有馬哲，石村貞夫， 東京図書，1987）、決定木（非特許文献１０：“AIによるデータ解析”，J.R.キンラン，トッパン，1995）、サポートベクターマシン（非特許文献１１：“サポートベクターマシン入門”，Nello Cristianini and John Shawe‐Taylor，共立出版，2005）等の分類的解析手法を用いることもできる。
いずれの場合も、例えば、説明変数を脳機能情報の時系列データρ(l,m,k,i)、外的基準をタスクとレストの時系列データとして、各種の判別関数、決定木等を用いてデータ分類する。本変形例では、その際に、ステップＳ４０において、データ評価量構成手段５は、各解析手法における通常の評価量に対して、ステップＳ３０で算出した結合度ベクトル
を組み込んだ評価量を構成する。そして、ステップＳ５０において、データ解析手段６は、各解析手法を用いてデータ解析を行う。上記のように評価量を変更すると、通常の解法が使えなくなる場合があるので、その際には、データ解析手段６は、最適なκを交差検証法を用いて決定し、GA等の確率的探索アルゴリズムによってデータ解析を行えばよい。
そして、ステップＳ６０において、画像生成手段７は、ステップＳ５０で得られた分類データを各種検定した際に有意であると判断されたボクセルの脳機能データρ(l,m,k)に基づく値を表示画像として生成し、ステップＳ７０において、画像表示手段８は、上記のようにして生成された画像を立体表示する。
これによって、第１の実施形態と同様の効果が得られる。
＜第１の実施形態の変形例５＞
第１の実施形態では、ステップＳ５０において、データ解析手段６は、式（２４）の線形回帰式を用いることによりノンパラメトリック回帰分析をしたわけだが、データ解析の手法として、ICA（Independent Component Analysis：独立成分分析）を用いることもできる（非特許文献１２：“詳解独立成分分析 信号解析の新しい世界”, Aapo Hyvarinen et al.，東京電機大学出版局 2005）。
本変形例では、ステップＳ４０において、データ評価量構成手段５は、ボクセル間結合度算出手段４で算出された結合度ベクトル
に基づき評価量を構成する。独立成分分析による脳機能情報の時系列データρ(l,m,k,i)の解析での「空間的な独立」という条件を、本変形例では「拡散テンソル情報にもとづく結合度を考慮した空間的な独立」というように変更する。独立成分分析の具体的手法はいくつかあるが、例えば、相互相関を最小化する手法の場合には、本変形例では、相互相関を最小化するのではなく、結合度に基づく評価量に近づけるというようにする。これにより、結合度の大きいボクセル間の相互相関は結合度に基づく評価量に近いような「独立度」になり、結合度の小さい（もしくは、結合度のない）ボクセル間は相互相関を最小にするというようになる。より具体的には、例えば、相関行列の非対角成分の２乗和を最小化するような評価関数中の２乗和の各項に、結合度に基づく評価量による重み係数をつける。この重み係数は、結合度が大きいほど小さくなるように構成する。
次に、ステップＳ５０において、データ解析手段６は、上記のようにして構成された評価量を考慮した独立成分分析を行う。その際、最適なκを交差検証法を用いて決定し、GA等の確率的探索アルゴリズムによってデータ解析を行う。
次に、ステップＳ６０において、画像生成手段７は、上記のようにして得られた独立成分を表示画像として生成する。ただし、上記の解析の結果得られる独立成分の数は多いので、全行ベクトルの中からタスク／レストと最も相関の高い行ベクトルのみを抽出する処理を行う。
そして、ステップＳ７０において、画像表示手段８は、上記のようにして生成された画像を立体表示する。
これによって、第１の実施形態と同様の効果が得られる。
なお、上記した独立成分分析以外にも、脳機能情報の時系列データから所定の主たる成分を抽出する手法として、主成分分析、因子分析、数量化III類等があるが、これらの解析手法についても、基本的には、上記と同様にして解析することができる。
[第２の実施形態]
図１０は、本発明の第２の実施形態に係る脳機能解析装置の概略構成を示したブロック図である。第２の実施形態の脳機能解析装置２０は、第１の実施形態の脳機能解析装置１０におけるボクセル間結合度算出手段４、データ評価量構成手段５、及びデータ解析手段６の代わりに、ボクセル間結合度算出手段２４、データ評価量構成手段２５、データ平滑化手段２００、及びデータ解析手段２６を設けた構成を成している。第１の実施形態と同一の構成には同一の符号を付し、以降においては、繰り返しの説明を省略する。なお、本実施形態においては、脳機能解析装置２０を、本実施形態の主要構成手段２０Ａと画像表示手段８及び記憶手段９とを一体型としたディスプレイコンソール形式としたが、画像表示手段８、或いは、画像表示手段８及び記憶手段９をそれぞれ独立した画像表示装置及び記憶装置として主要構成手段２０Ａから切り離した構成としてもよい。
ボクセル間結合度算出手段２４は、拡散テンソルデータ取得手段２で取得された拡散テンソルデータD(l,m,k)（式（９））から第１の実施形態に詳述したようにして平均拡散度ベクトル
（式（１２））を算出した後、該平均拡散度ベクトル
の各成分D’l(l,m,k)，D’m(l,m,k)，D’k(l,m,k)の大きさに応じてデータ前処理手段３で前処理された脳機能情報の時系列データρ(l,m,k,i)を平滑化する際の隣接ボクセル間の平滑度を決める結合度ベクトル
を算出する。
データ平滑化手段２００は、ボクセル間結合度算出手段２４で算出された結合度ベクトル
を用いて、データ前処理手段３で前処理された脳機能情報の時系列データρ(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)から第１の実施形態で詳述したようにして平均拡散度ベクトル
を算出する。次いで、該平均拡散度ベクトル
の各成分D’l(l,m,k)，D’m(l,m,k)，D’k(l,m,k)の大きさに応じてデータ前処理手段３で前処理された脳機能情報の時系列データρ(l,m,k,i)を平滑化する際の隣接ボクセル間の平滑度を決める結合度ベクトル
を、例えば、
として算出する。
次に、ステップＳ１３０において、データ平滑化手段２６は、ステップＳ１２０において算出された結合度ベクトル
を用いて、ステップＳ１１０において前処理された脳機能情報の時系列データρ(l,m,k,i)を、例えば、
のように異方的に平滑化する。ここで、w及びWは重み因子であり、
を満たすように規格化する。
図１２は、本実施形態に係るデータ平滑化の手法を例示したものであり、ある２次元スライス画像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を満たすような状況を想定する。このとき、結合度ベクトル
は、式（３１）から
となり、位置(l,m,k)のボクセル上の脳機能情報の時系列データρ(l,m,k,i)=5は、式（３６）から、
のように平滑化されることがわかる。この処理を該２次元スライス画像Skの全ボクセル上の脳機能情報の時系列データρ(l,m,k,i)に施すことによって、神経繊維の走行方向に沿った方向（つまり、平均拡散度ベクトル
の成分D’l(l,m,k)，D’m(l,m,k)，D’k(l,m,k)のうち、最も大きい成分の方向：この例では、X方向）からの寄与が大きいように平滑化することができる。
また、これ以外にも、式（５）の箇所で説明したように、最大固有値λMに対応する固有ベクトル
の向きと合致する２個のボクセルとのみ平滑化してもよい。
次に、ステップＳ１４０において、データ解析手段２６は、ステップＳ１４０において平滑化された脳機能情報の時系列データ
に対して、SPM等の解析手法を用いてデータ解析を施す。
次に、ステップＳ１５０において、画像生成手段７は、データ解析手段２６によって解析された結果をグレイスケールやフルカラー画像等の色調画像を生成する。
そして、ステップＳ１６０において、画像表示手段８は、上記のようにして生成された画像を立体表示する。
このように、第２の実施形態においては、MRI装置から取得した脳機能情報の時系列データを、同じくMRI装置から取得した拡散テンソルデータから算出された結合度ベクトルに基づいて平滑化処理したことによって、脳の活性部位間の結合構造をも考慮した脳の活性部位の特定が可能となる。
なお、式（３１）の結合度ベクトル
や式（３２）の平滑化処理のアルゴリズムはあくまでも一例であって、その他の実施例も可能である。さらに、本実施形態におけるデータ平滑化処理の手法は、特定のデータ解析手法に依存しないことは上記の説明から明らかである。
＜第２の実施形態の変形例＞
第２の実施形態においては、ステップＳ１３０において脳機能情報の時系列データρ(l,m,k,i)に対して拡散テンソルデータD(l,m,k)から算出された結合度ベクトル
を用いて平滑化処理を施した後に、ステップＳ１４０においてSPM等の解析手法を用いてデータ解析を行ったが、本発明の目的を達成するにあたって、これらのステップを逆にするような手法も考えられる。そこで、次に、上記した第２の実施形態の変形例について説明する。なお、重複を避けるため、以降においては第２の実施形態と異なる部分についてのみ説明を加える。
図１３は、第２の実施形態の変形例に係る脳機能解析装置の概略構成を示したブロック図である。本変形例に係る脳機能解析装置２０’の構成は、第２の実施形態の構成と基本的には同じであるが、データ平滑化手段２００とデータ解析手段２６とが入れ替わった構成になっている。これによって、脳機能解析装置２０’によって行われるデータ解析の処理手順が第２の実施形態と異なる。なお、本実施形態においては、脳機能解析装置２０’を、本実施形態の主要構成手段２０Ａ’と画像表示手段８及び記憶手段９とを一体型としたディスプレイコンソール形式としたが、画像表示手段８、或いは、画像表示手段８及び記憶手段９をそれぞれ独立した画像表示装置及び記憶装置として主要構成手段２０Ａ’から切り離した構成としてもよい。
図１４は、図１３に示した脳機能解析装置２０’によって行われるデータ解析の処理手順の一例を示したフローチャートである。本変形例では、ステップＳ１２１において、まず、データ解析手段２６’は、データ前処理手段３によって前処理された脳機能情報の時系列データρ(l,m,k,i)をSPM等の解析手法を用いてデータ解析する。次に、ステップＳ１３１において、ボクセル間結合度算出手段２４は、第２の実施形態と同様にして、式（３１）で示したような結合度ベクトル
を算出する。次いで、ステップＳ１４１において、データ平滑化手段２００’は、データ解析手段２６’によって解析された結果である脳機能データρ(l,m,k)に対して、第２の実施形態と同様にして、式（３２）を用いて平滑化処理を行う。次いで、ステップＳ１５０において、画像生成手段７は、上記のようにして得られた平滑化処理後のデータ
に対して、所定の処理規則に従って、グレイスケールやフルカラー画像等の色調画像を生成する。そして、ステップＳ１６０において、画像表示手段８は、上記のようにして生成された画像を立体表示する。
このようにして、本変形例においても、第２の実施形態と同様の効果が得られる。
[第３の実施形態]
図１５は、本発明の第３の実施形態に係る脳機能解析装置の概略構成を示したブロック図である。第３の実施形態の脳機能解析装置３０は、第１の実施形態の脳機能解析装置１０におけるボクセル間結合度算出手段４、データ評価量構成手段５、及びデータ解析手段６の代わりに、ボクセル間結合度算出手段３４、クラスタリング手段３００、及びデータ解析手段３６を設けた構成を成している。第１の実施形態と同一の構成には同一の符号を付し、以降においては、繰り返しの説明を省略する。なお、本実施形態においては、脳機能解析装置３０を、本実施形態の主要構成手段３０Ａと画像表示手段８及び記憶手段９とを一体型としたディスプレイコンソール形式としたが、画像表示手段８、或いは、画像表示手段８及び記憶手段９をそれぞれ独立した画像表示装置及び記憶装置として主要構成手段３０Ａから切り離した構成としてもよい。
ボクセル間結合度算出手段３４は、拡散テンソルデータ取得手段２で取得された拡散テンソルデータD(l,m,k)（式（９））から第１の実施形態で詳述したようにして平均拡散度ベクトル
（式（１２））を算出した後、該平均拡散度ベクトル
の各成分D’l(l,m,k)，D’m(l,m,k)，D’k(l,m,k)の大きさに応じてボクセル同士をクラスタリングするための指標として結合度ベクトル
を算出する。
クラスタリング手段３００は、ボクセル間結合度算出手段３４によって算出された結合度ベクトル
を用いてボクセルをクラスタリングする。
データ解析手段３６は、クラスタリング手段３００でクラスタリングされたボクセル集団（クラスター）を単位として、データ前処理手段３で前処理された脳機能情報の時系列データρ(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)から第１の実施形態で詳述したようにして平均拡散度ベクトル
を算出する。次いで、該平均拡散度ベクトル
の各成分D’l(l,m,k)，D’m(l,m,k)，D’k(l,m,k)の大きさに応じてボクセル同士をクラスタリングする際の指標となる結合度ベクトル
を、例えば、
として算出する。ここで、
はクラスタリングの閾値を表す定数ベクトルであり、
と書ける。
次に、ステップＳ２３０において、クラスタリング手段３００は、ステップＳ２２０において算出された結合度ベクトル
の各成分Sl(l,m,k)，Sm(l,m,k),Sk(l,m,k)の値が正になるボクセルを、それぞれX方向、Y方向、Z方向に対してクラスタリングする。図１７は、図１５に示したクラスタリング手段３００によって行われるクラスタリングの手法を例示した図であり、ある２次元スライス画像Skの１６個のボクセルを表している。図１７（Ａ）における隣接ボクセル間の矢印の値は、X方向又はY方向の平均拡散度ベクトルの値を表している。式（３７）の定数ベクトル
の各成分の値をそれぞれ、例えば、hl=hm=3，hk=0とすると、結合度ベクトル
の成分Sl(l,m,k)，Sm(l,m,k)の各値は、図１７（Ｂ）のようになる。このとき、X方向、Y方向の各方向において、結合度ベクトル
の成分Sl(l,m,k)，Sm(l,m,k)が共に正となる隣接ボクセルを順次結んでいくことによって、最終的に、点線で囲んだボクセルが一つのクラスターとして形成されることになる。結合度ベクトル
が平均拡散度ベクトルD’(l,m,k)の関数となっていることから、このようなクラスタリングは神経線維の走行方向を反映したものとなっている。
次に、ステップＳ２４０において、データ解析手段３６は、ステップＳ２３０においてクラスタリングされたボクセル集団（クラスター）毎に、ステップＳ２１０において前処理された脳機能情報の時系列データρ(l,m,k,i)を、SPM等の解析方法を用いてデータ解析する。その際には、例えば、一つのクラスター内での脳機能情報の時系列データを平均化する等の処理が必要となる。
次に、ステップＳ２５０において、画像生成手段７は、ステップＳ２４０においてデータ解析された後のデータ
に対して、所定の処理規則に従って、グレイスケールやフルカラー画像等の色調画像を生成する。
そして、ステップＳ２６０において、画像表示手段８は、ステップＳ２５０において生成された画像を立体表示する。
このように、第３の実施形態においては、MRI装置から取得した拡散テンソルデータに基づいてボクセルをクラスタリングした後に、同じくMRI装置から取得した脳機能情報の時系列データをクラスター毎に解析したことによって、脳の活性部位間の結合構造をも考慮した脳の活性部位の特定が可能となる。
なお、式（３６）の結合度ベクトル
はあくまでも一例であってその他の実施例も可能である。また、式（５）の箇所で説明したように、最大固有値λMに対応する固有ベクトル
の向きと合致する２個のボクセルとのみクラスタリングしてもよい。さらに、本実施形態におけるクラスタリング処理の手法は、特定のデータ解析手法に依存しないことは上記の説明から明らかである。
＜第３の実施形態の変形例１＞
第３の実施形態においては、ステップＳ２３０において、拡散テンソルデータD(l,m,k)から算出された結合度ベクトル
を用いてボクセルをクラスタリングした後に、ステップＳ２４０においてクラスター毎に脳機能情報の時系列データρ(l,m,k,i)をSPM等の解析手法を用いてデータ解析したが、本発明の目的を達成するにあたって、これらのステップを逆にするような手法も考えられる。そこで、次に、上記した第３の実施形態の変形例について説明する。なお、重複を避けるため、以降においては第３の実施形態と異なる部分についてのみ説明を加える。
図１８は、第３の実施形態の変形例に係る脳機能解析装置の概略構成を示したブロック図である。本変形例に係る脳機能解析装置３０’の構成は、第３の実施形態の構成と基本的には同じであるが、クラスタリング手段３００とデータ解析手段３６とが入れ替わった構成になっている。これによって、脳機能解析装置３０’によって行われるデータ解析の処理手順が第３の実施形態と異なる。なお、本実施形態においては、脳機能解析装置３０’を、本実施形態の主要構成手段３０Ａ’と画像表示手段８及び記憶手段９とを一体型としたディスプレイコンソール形式としたが、画像表示手段８、或いは、画像表示手段８及び記憶手段９をそれぞれ独立した画像表示装置及び記憶装置として主要構成手段３０Ａ’から切り離した構成としてもよい。
図１９は、図１８に示した脳機能解析装置３０’によって行われるデータ解析の処理手順の一例を示したフローチャートである。本変形例では、ステップＳ２２１において、まず、データ解析手段３６’は、データ前処理手段３によって前処理された脳機能情報の時系列データρ(l,m,k,i)をSPM等の解析手法でデータ解析する。次に、ステップＳ２３１において、ボクセル間結合度算出手段３４は、第３の実施形態と同様にして、式（３６）で示したような結合度ベクトル
を算出する。次いで、ステップＳ２４１において、クラスタリング手段３００’は、データ解析手段３６’によって解析された脳機能データρ(l,m,k)に対して、第３の実施形態と同様にして、図１７に示したようなボクセルのクラスタリングを行う。次いで、ステップＳ２５０において、画像生成手段７は、上記のようにしてクラスタリングされた脳機能データ
に対して、所定の処理規則に従って、グレイスケールやフルカラー画像等の色調画像を生成する。そして、ステップＳ２６０において、画像表示手段８は、上記のようにして生成された画像を立体表示する。
このようにして、本変形例においても、第３の実施形態と同様の効果が得られる。
＜第３の実施形態の変形例２＞
第３の実施形態の変形例１においては、データ前処理手段３によって前処理された脳機能情報の時系列データρ(l,m,k,i)をSPM等の解析手法でデータ解析した後に、ボクセル間結合度算出手段３４で算出された結合度ベクトル
を用いてデータ解析後の脳機能データρ(l,m,k)をクラスタリングしたが、データ解析の手法として分類的解析手法を用いることも考えられる。そこで、本変形例では、その一例として、データ解析の手法として分類的解析手法の一つである決定木を用いた場合について説明する。例えば、図４に示したように、説明変数（属性）を脳機能情報の時系列データρ(l,m,k)、外的基準（クラス）をタスクとレストの時系列データとして決定木の手法を用いた場合、不要な属性（ボクセル）は除去され、必要な属性（ボクセル）のみが残る。このようにして残ったボクセルを結合度ベクトル
を用いてクラスタリングする。例えば、（１）あるボクセルを選び、（２）そのボクセルの周囲で結合度が一定値以上のボクセルをクラスタリングし、（３）（２）の処理を、周囲に結合度が一定値以上のボクセルがなくなるまで続け、（４）ボクセルが残っていれば、あるボクセルを選び、（２）の処理にもどる。
次いで、ステップＳ２５０において、画像生成手段７は、上記のようにしてクラスタリングされた脳機能データ
に対して、所定の処理規則に従って、グレイスケールやフルカラー画像等の色調画像を生成する。そして、ステップＳ２６０において、画像表示手段８は、上記のようにして生成された画像を立体表示する。
このようにして、本変形例においても、第３の実施形態と同様の効果が得られる。
[第４の実施形態]
図２０は、第４の実施形態に係る脳機能解析装置の概略構成を示したブロック図である。第４の実施形態の脳機能解析装置４０は、第１の実施形態の脳機能解析装置１０におけるボクセル間結合度算出手段４、データ評価量構成手段５、及びデータ解析手段６の代わりに、ボクセル間結合度算出手段４４、データ検定手段４００を設けた構成を成している。第１の実施形態と同一の構成には同一の符号を付し、以降においては、繰り返しの説明を省略する。なお、本実施形態においては、脳機能解析装置４０を、本実施形態の主要構成手段４０Ａと画像表示手段８及び記憶手段９とを一体型としたディスプレイコンソール形式としたが、画像表示手段８、或いは、画像表示手段８及び記憶手段９をそれぞれ独立した画像表示装置及び記憶装置として主要構成手段４０Ａから切り離した構成としてもよい。
ボクセル間結合度算出手段４４は、データ前処理手段３で前処理された脳機能情報の時系列データρ(l,m,k,i)をボクセル毎に各種データ検定（例えば、ｔ検定、ｆ検定、又はｚ検定）する際の基準値（検定値や有意水準など）を調整するための結合度ベクトル
を、式（１２）の平均拡散度ベクトル
の関数として、拡散テンソルデータ取得手段２で取得された拡散テンソルデータD(l,m,k)（式（９））から算出する。
データ検定手段４００は、ボクセル間結合度算出手段４４で算出された結合度ベクトル
によって該基準値を修正し、修正した基準値に基づき、前処理後の脳機能データρ(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)をボクセル毎に、例えば、ｚ検定する際の基準値（検定値や有意水準など）を調整するための結合度ベクトル
を、式（１２）の平均拡散度ベクトル
の関数として、ステップＳ３１０で取得した拡散テンソルデータD(l,m,k)から算出する。
次に、ステップＳ３３０において、データ検定手段４００は、ステップＳ３１０において前処理された脳機能情報の時系列データρ(l,m,k,i)をボクセル毎にｚ検定を行いｚスコアを算出する。そして、同ステップにおいて、算出したｚスコアが大きなボクセル（例えば、ボクセル(l,m,k)）を基点にして、その周辺（２６方位）のボクセル（例えば、ボクセル(l+1,m,k)）における有意水準の値a(l+1,m,k)を、ステップＳ３２０で算出した結合度ベクトル
に応じて、例えば、
のように下げる。ここで、αは重み係数である。
次に、ステップＳ３４０において、画像生成手段７は、ステップＳ３３０において修正された有意水準に基づき、有意であると判断されたボクセル上の脳機能データρ(l,m,k)に基づく値を画像データとして生成する。
そして、ステップＳ３５０において、画像表示手段８は、ステップＳ３４０において生成された画像データを立体表示する。
なお、上記においては、拡散テンソルデータ取得手段２で取得した拡散テンソルデータD(l,m,k)から算出した式（３８）の結合度ベクトル
を用いて検定の有意水準の値を調整したが、その代わりに、検定値を調整してもよい。
例えば、タスク画像とレスト画像との間でボクセル毎にz検定を行う場合を考えてみる。タスク画像とは、あるタスクを行っているときの画像（例えば、図２の太線で表された１〜３のスライス画像）のことをいい、レスト画像とは、該タスクを行っていないときの画像（例えば、図２の太線で表された５〜７のスライス画像）のことをいう。このとき、ボクセル(l,m,k)におけるzスコア（z検定の検定値）は、通常、
と表される。ここで、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スコアの値が大きいボクセルがタスクに関連している活性部位であると考えることができる。
そこで、式（４０）を
のように変更して、ステップＳ３３０において、b(l,m,k)を、ステップＳ３２０で算出した式（３８）の結合度ベクトル
を用いて、例えば、
と修正する。ここで、βは重み係数である。
次いで、ステップＳ３４０において、画像生成手段７は、全ボクセルにおいて共通の有意水準に基づき、有意であると判断されたボクセル上の脳機能データρ(l,m,k)に基づく値を画像データとして生成する。そして、ステップＳ３５０において、画像表示手段８は、ステップＳ３４０において生成された画像データを立体表示する。
このように、第４の実施形態においては、MRI装置から取得した拡散テンソルデータに基づいて、同じくMRI装置から取得した脳機能情報の時系列データを検定する際の基準値（検定値又は有意水準）の値をボクセル毎に調整したことによって、脳の活性部位間の結合構造をも考慮した脳の活性部位の特定が可能となる。
＜第４の実施形態の変形例＞
第４の実施形態では、脳機能情報の時系列データρ(l,m,k,i)をボクセル毎に各種データ検定（例えば、Ｔ検定、Ｆ検定、又はＺ検定）する際に、拡散テンソルデータ取得手段２で取得した拡散テンソルデータD(l,m,k)から算出した結合度ベクトル
を用いて該データ検定の基準値（検定値や有意水準など）を調整したわけだが、本発明の目的を達成するにあたって、脳機能情報の時系列データρ(l,m,k,i)をボクセル毎に各種の解析手法でデータ解析したときに、該解析手法が統計的に有意であるか否かをデータ検定（例えば、Ｔ検定、Ｆ検定、又はＺ検定）するような状況も考えられる。そこで、本変形例として、その際の基準値（検定値や有意水準など）を、結合度ベクトル
を用いて調整することを考える。なお、重複を避けるため、以降においては第４の実施形態と異なる部分についてのみ説明を加える。
図２２は、第４の実施形態の変形例に係る脳機能解析装置の概略構成を示した概念図である。本変形例に係る脳機能解析装置４０’の構成は、第４の実施形態の構成のデータ検定手段４００の代わりに、データ検定手段４００’とデータ解析手段４６’とを加えた構成を成す。
データ解析手段４６’は、脳機能データ取得手段１で取得され、データ前処理手段３で前処理された脳機能情報の時系列データρ(l,m,k,i)を、従来の分類、回帰、相関等の解析手法を用いてボクセル毎にデータ解析する。
データ検定手段４００’は、データ解析手段４６’で行われた解析手法が統計的に有意であるか否かを各種データ検定（例えば、Ｔ検定、Ｆ検定、又はＺ検定）に基づき検定する際の基準値（検定値や有意水準など）を、結合度ベクトル
を用いて調整する。
図２３は、図２２に示した脳機能解析装置によって行われるデータ解析の処理手順の一例を示したフローチャ−トである。
ステップＳ３６０において、データ解析手段４６’は、脳機能データ取得手段１で取得され、データ前処理手段３で前処理された脳機能情報の時系列データρ(l,m,k,i)を、従来の分類、回帰、相関等の解析手法を用いてボクセル毎にデータ解析する。
１．分類的解析手法（判別分析、決定木、サポートベクターマシン等）
１−１：判別分析の場合
説明変数を脳機能情報の時系列データρ(l,m,k,i)、外的基準をタスクとレストの時系列データとして、各種の判別関数を用いてデータ分類する。
１−２：決定木の場合
説明変数（属性）をタスクとレストの時系列データ、外的基準として脳機能情報の時系列データρ(l,m,k,i)を離散化して幾つかのクラスを作り、決定木を用いてデータ分類する。
２．回帰的解析手法（各種の関数（線形、非線形）を用いた回帰分析）
２−１：線形回帰分析の場合
説明変数をタスクとレスト、被説明変数を脳機能情報の時系列データρ(l,m,k,i)として、線形回帰式を用いてデータ解析する。
３．相関的解析手法（推定モデル（関数）と実際のデータとの相関を取る手法）
予め設定した、モデル（例えば、刺激を入力とし、脳機能情報の時系列データを出力とするような関数）と実際の脳機能情報の時系列データρ(l,m,k,i)との相関を取ることによって、データ解析する。
ついで、ステップＳ３２０（第４の実施形態と同様）を経て、ステップＳ３３１において、データ検定手段４００’は、ステップＳ３６０において得られた結果が統計的に有意であるか否かを検定する。このとき、検定の基準値（検定値や有意水準など）を結合度ベクトル
（式（３８））を用いて調整する。例えば、分類的解析手法で得られた結果に対してf検定を行う場合には、ステップＳ３６０において得られた結果に対してボクセル毎にf検定を行いf値を算出する。そして、同ステップにおいて、算出したf値が大きなボクセル（例えば、ボクセル(l,m,k)）を基点にして、その周辺（２６方位）のボクセル（例えば、ボクセル(l+1,m,k)）における有意水準の値a(l+1,m,k)を、ステップＳ３２０で算出した結合度ベクトル
に応じて、例えば、式（３９）を用いて調整する。回帰的解析手法の場合には、閾値（有意水準に相当）の値を上記のようにして調整すればよいし、相関的解析手法の場合には、相関の閾値を上記のようにして調整すればよい。
次に、ステップＳ３４０において、画像生成手段７は、ステップＳ３３０において修正された有意水準に基づき、有意であると判断されたボクセル上の脳機能データρ(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とトラクトグラフィーとを組み合わせた解析手法においては、神経繊維（この場合、弓状束）の検出にあたり、解析者が神経繊維の追跡開始点と終了点とを予め指定しなければならない。したがって、神経線維の追跡開始点と終了点とが予め決定できないような場合（例えば、灰白質の活性領域が予めわかっていないようなタスク）にはこの手法は有効でないことになる。それに対し、本発明の第１の実施形態に係る脳機能解析方法によれば、神経線維の追跡開始点と終了点とが予め指定できない場合に対しても、図２４Ｄに示したように、神経線維の活性化の検出が可能である。
このように、本発明によれば、脳の活性部位を特定できるのはもちろんのこと、解析者が予め神経線維の開始点及び終了点を指定しなくとも、神経線維（活性部位間の結合構造）を特定することが可能となり、さらには神経繊維の活性化（神経線維におけるパルス伝達）をも検出可能となるので、認知症、失語症、統合失調症等の高次脳機能障害の診断及び治療に対して極めて有効な脳機能解析方法、及び脳機能解析プログラムを提供することができる。

Claims (27)

  1. ボクセル単位で取得された、脳内のプロトンの拡散度を特定可能とする拡散テンソルデータに基づき隣接ボクセル間の結合度の評価量を構成し、
    前記隣接ボクセル間の結合度の評価量に基づき、前記ボクセル単位で取得された、脳の活性部位を特定可能とする脳機能データの解析を行う、
    ことを含む、脳機能データ解析方法。
  2. 前記脳機能データの解析を行うことは、前記脳機能データ及びタスクの一方を説明変数とし、前記脳機能データ及び前記タスクの他方を被説明変数とする回帰分析により、前記脳機能データのボクセル毎の評価量に前記隣接ボクセル間の結合度の評価量を組み込んだ評価量の最小値及び最大値の一方を算出することを含む、請求項１記載の脳機能データ解析方法。
  3. 前記脳機能データの解析を行うことは、前記脳機能データの検定の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを含む、請求項１記載の脳機能データ解析方法。
  4. 前記脳機能データの解析を行うことは、前記脳機能データの分類の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを含む、請求項１記載の脳機能データ解析方法。
  5. 前記脳機能データの解析を行うことは、前記脳機能データと所定のモデルとの相関の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを含む、請求項１記載の脳機能データ解析方法。
  6. 前記脳機能データの解析を行うことは、前記脳機能データから主要な成分を抽出する抽出の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを含む、請求項１記載の脳機能データ解析方法。
  7. 前記脳機能データの解析を行うことは、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データを平滑化することを含む、請求項１記載の脳機能データ解析方法。
  8. 前記脳機能データの解析を行うことは、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データをクラスタリングすることを含む、請求項１記載の脳機能データ解析方法。
  9. 前記拡散テンソルデータに基づき、前記脳機能データから特定される活性部位ボクセルの脳機能データの値を他のボクセルの脳機能データの値に伝播させる前処理を施すことをさらに含む、請求項１記載の脳機能データ解析方法。
  10. コンピュータを、
    脳の活性部位を特定可能とする脳機能データをボクセル単位で取得する脳機能データ取得手段と、
    前記脳内のプロトンの拡散度を特定可能とする拡散テンソルデータをボクセル単位で取得する拡散テンソルデータ取得手段と、
    前記拡散テンソルデータに基づき隣接ボクセル間の結合度の評価量を構成するデータ評価量構成手段と、
    前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データの解析を行うデータ解析手段と、
    して機能させることを特徴とする脳機能解析プログラム。
  11. 請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記脳機能データ及びタスクの一方を説明変数とし、前記脳機能データ及び前記タスクの他方を被説明変数とする回帰分析する手段であって、前記脳機能データのボクセル毎の評価量に前記隣接ボクセル間の結合度の評価量を組み込んだ評価量の最小値及び最大値の一方を算出することを特徴とする脳機能解析プログラム。
  12. 請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記脳機能データを検定する手段であって、前記検定の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析プログラム。
  13. 請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記脳機能データを分類する手段であって、前記分類の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析プログラム。
  14. 請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記脳機能データと所定のモデルとの相関をとる手段であって、前記相関の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析プログラム。
  15. 請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記脳機能データから主要な成分を抽出する手段であって、前記抽出の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析プログラム。
  16. 請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析は、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データを平滑化することを特徴とする脳機能解析プログラム。
  17. 請求項１０に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データをクラスタリングすることを特徴とする脳機能解析プログラム。
  18. 請求項１０乃至１７のいずれかに記載の脳機能解析プログラムにおいて、前記コンピュータを、前記拡散テンソルデータ取得手段によって取得した拡散テンソルデータに基づき、前記脳機能データ取得手段によって取得した脳機能データから特定される活性部位ボクセルの脳機能データの値を他のボクセルの脳機能データの値に伝播させる前処理を施すデータ前処理手段として機能させることを特徴とする脳機能解析プログラム。
  19. 脳の活性部位を特定可能とする脳機能データをボクセル単位で取得する脳機能データ取得手段と、
    前記脳内のプロトンの拡散度を特定可能とする拡散テンソルデータをボクセル単位で取得する拡散テンソルデータ取得手段と、
    前記拡散テンソルデータに基づき隣接ボクセル間の結合度の評価量を構成するデータ評価量構成手段と、
    前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データの解析を行うデータ解析手段と、
    を備えたことを特徴とする脳機能解析装置。
  20. 請求項１９に記載の脳機能解析装置において、前記データ解析手段は、前記脳機能データ及びタスクの一方を説明変数とし、前記脳機能データ及び前記タスクの他方を被説明変数とする回帰分析する手段であって、前記脳機能データのボクセル毎の評価量に前記隣接ボクセル間の結合度の評価量を組み込んだ評価量の最小値及び最大値の一方を算出することを特徴とする脳機能解析装置。
  21. 請求項１９に記載の脳機能解析装置において、前記データ解析手段は、前記脳機能データを検定する手段であって、前記検定の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析装置。
  22. 請求項１９に記載の脳機能解析装置において、前記データ解析手段は、前記脳機能データを分類する手段であって、前記分類の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析装置。
  23. 請求項１９に記載の脳機能解析装置において、前記データ解析手段は、前記脳機能データと所定のモデルとの相関をとる手段であって、前記相関の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析装置。
  24. 請求項１９に記載の脳機能解析装置において、前記データ解析手段は、前記脳機能データから主要な成分を抽出する手段であって、前記抽出の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析装置。
  25. 請求項１９に記載の脳機能解析装置において、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データを平滑化するデータ平滑化手段をさらに備えたことを特徴とする脳機能解析装置。
  26. 請求項１９に記載の脳機能解析装置において、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データをクラスタリングするデータクラスタリング手段をさらに備えたことを特徴とする脳機能解析装置。
  27. 請求項１９乃至２６のいずれかに記載の脳機能解析装置において、前記拡散テンソルデータ取得手段によって取得した拡散テンソルデータに基づき、前記脳機能データ取得手段によって取得した脳機能データから特定される活性部位ボクセルの脳機能データの値を他のボクセルの脳機能データの値に伝播させる前処理を施すデータ前処理手段をさらに備えたことを特徴とする脳機能解析装置。
JP2007539917A 2005-10-12 2006-10-06 脳機能データ解析方法、脳機能解析装置及び脳機能解析プログラム Active JP4308871B2 (ja)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
JP2005298236 2005-10-12
JP2005298236 2005-10-12
PCT/JP2006/320078 WO2007043462A1 (ja) 2005-10-12 2006-10-06 脳機能解析方法および脳機能解析プログラム

Publications (2)

Publication Number Publication Date
JPWO2007043462A1 JPWO2007043462A1 (ja) 2009-04-16
JP4308871B2 true JP4308871B2 (ja) 2009-08-05

Family

ID=37942705

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2007539917A Active JP4308871B2 (ja) 2005-10-12 2006-10-06 脳機能データ解析方法、脳機能解析装置及び脳機能解析プログラム

Country Status (6)

Country Link
US (1) US8059879B2 (ja)
EP (1) EP1946701B1 (ja)
JP (1) JP4308871B2 (ja)
KR (1) KR100971936B1 (ja)
CN (1) CN101287410B (ja)
WO (1) WO2007043462A1 (ja)

Families Citing this family (33)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2009146388A1 (en) 2008-05-28 2009-12-03 The Trustees Of Columbia University In The City Of New York Voxel-based methods for assessing subjects using positron emission tomography
WO2010048438A1 (en) * 2008-10-22 2010-04-29 The Trustees Of Columbia University In The City Of New York Images of language-sensitive neurocircuitry as a diagnostic for autism
US8666475B2 (en) 2008-10-22 2014-03-04 The Trustees Of Columbia University In The City Of New York Images of language-sensitive neurocircuitry as a diagnostic for autism
CN101912263B (zh) * 2010-09-14 2011-12-14 北京师范大学 基于脑功能网络成分检测的实时功能磁共振数据处理系统
US8467585B2 (en) * 2010-10-21 2013-06-18 Genenal Electric Company Methods and apparatus to analyze computed tomography scan data
US20120163689A1 (en) * 2010-12-28 2012-06-28 Boettger Joachim Method and device for visualizing human or animal brain segments
CN102973279B (zh) * 2012-12-18 2014-09-17 哈尔滨工业大学 独立成分分析联合最小二乘法的近红外脑机接口的信号检测方法
KR101401969B1 (ko) * 2013-02-20 2014-06-30 한양대학교 산학협력단 Ｍｒｉ 시스템을 이용하여 대상체 내의 신경 섬유에 관한 섬유 구조 정보를 획득하기 위한 방법 및 장치
CN104207776A (zh) * 2014-08-22 2014-12-17 南昌大学 一种综合性磁共振成像装置及方法
JP6163471B2 (ja) * 2014-09-17 2017-07-12 株式会社日立製作所 磁気共鳴イメージング装置
CN104323776B (zh) * 2014-10-23 2016-11-16 中国科学院深圳先进技术研究院 脑功能磁共振成像方法和系统
CN104571533B (zh) * 2015-02-10 2018-03-13 北京理工大学 一种基于脑机接口技术的装置和方法
CN105069307B (zh) * 2015-08-19 2018-04-10 大连理工大学 一种结合独立成分分析与移不变规范多元分解的多被试功能核磁共振成像数据分析方法
JP6707330B2 (ja) * 2015-09-10 2020-06-10 キヤノンメディカルシステムズ株式会社 画像処理装置及び磁気共鳴イメージング装置
CN105816170B (zh) * 2016-05-10 2019-03-01 广东省医疗器械研究所 基于可穿戴式nirs-eeg的精神分裂症早期检测评估系统
CN106940803B (zh) * 2017-02-17 2018-04-17 平安科技（深圳）有限公司 相关变量识别方法和装置
TWI651688B (zh) * 2017-03-17 2019-02-21 長庚大學 利用磁振造影影像預測神經疾病的臨床嚴重度的方法
TWI639830B (zh) * 2017-03-17 2018-11-01 長庚大學 Method for identifying neurological diseases using magnetic resonance imaging images
WO2019060298A1 (en) 2017-09-19 2019-03-28 Neuroenhancement Lab, LLC METHOD AND APPARATUS FOR NEURO-ACTIVATION
US11717686B2 (en) 2017-12-04 2023-08-08 Neuroenhancement Lab, LLC Method and apparatus for neuroenhancement to facilitate learning and performance
US11478603B2 (en) 2017-12-31 2022-10-25 Neuroenhancement Lab, LLC Method and apparatus for neuroenhancement to enhance emotional response
US11364361B2 (en) 2018-04-20 2022-06-21 Neuroenhancement Lab, LLC System and method for inducing sleep by transplanting mental states
EP3581090A1 (en) * 2018-06-11 2019-12-18 Koninklijke Philips N.V. Electrical properties tomography mapping of conductivity changes
WO2020056418A1 (en) 2018-09-14 2020-03-19 Neuroenhancement Lab, LLC System and method of improving sleep
EP3657393A1 (en) * 2018-11-20 2020-05-27 Koninklijke Philips N.V. Determination of a further processing location in magnetic resonance imaging
JP7321703B2 (ja) * 2018-12-19 2023-08-07 富士フイルムヘルスケア株式会社 画像処理装置、及び、磁気共鳴イメージング装置
CN109830286B (zh) * 2019-02-13 2022-09-30 四川大学 基于非参数统计的脑功能磁共振编码能量成像方法
US11786694B2 (en) 2019-05-24 2023-10-17 NeuroLight, Inc. Device, method, and app for facilitating sleep
CN111081351B (zh) * 2019-12-02 2024-02-20 北京优脑银河科技有限公司 脑功能图谱的绘制方法和系统
US20210170180A1 (en) 2019-12-09 2021-06-10 Washington University Systems and methods of precision functional mapping-guided personalized neuromodulation
US11151456B1 (en) 2021-01-08 2021-10-19 Omniscient Neurotechnology Pty Limited Predicting brain data using machine learning models
US11666266B2 (en) * 2021-10-05 2023-06-06 Omniscient Neurotechnology Pty Limited Source localization of EEG signals
CN116523857B (zh) * 2023-04-21 2023-12-26 首都医科大学附属北京友谊医院 基于弥散张量图像的听力状态预测装置、方法

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH0947438A (ja) * 1995-08-07 1997-02-18 Hitachi Ltd 活性化領域同定法
JP2004081657A (ja) 2002-08-28 2004-03-18 Ge Medical Systems Global Technology Co Llc 繊維状画像抽出方法、画像処理装置および磁気共鳴撮像システム
US7627370B2 (en) * 2004-10-20 2009-12-01 Marks Donald H Brain function decoding process and system
US7894903B2 (en) * 2005-03-24 2011-02-22 Michael Sasha John Systems and methods for treating disorders of the central nervous system by modulation of brain networks

Also Published As

Publication number Publication date
CN101287410A (zh) 2008-10-15
EP1946701A1 (en) 2008-07-23
WO2007043462A1 (ja) 2007-04-19
EP1946701B1 (en) 2013-08-21
US8059879B2 (en) 2011-11-15
JPWO2007043462A1 (ja) 2009-04-16
EP1946701A4 (en) 2009-11-18
KR100971936B1 (ko) 2010-07-23
KR20080058475A (ko) 2008-06-25
US20090279762A1 (en) 2009-11-12
CN101287410B (zh) 2013-09-11

Similar Documents

Publication Publication Date Title
JP4308871B2 (ja) 脳機能データ解析方法、脳機能解析装置及び脳機能解析プログラム
Tohka et al. Automatic independent component labeling for artifact removal in fMRI
Sui et al. An ICA-based method for the identification of optimal FMRI features and components using combined group-discriminative techniques
US9271679B2 (en) Method and apparatus for processing medical image signal
Loued-Khenissi et al. An overview of functional magnetic resonance imaging techniques for organizational research
US20170238879A1 (en) Method of Analyzing the Brain Activity of a Subject
EP2147330B1 (en) Image processing method
US10736538B2 (en) Method and computer differentiating correlation patterns in functional magnetic resonance imaging
Luckett et al. Resting state network mapping in individuals using deep learning
Turner Where matters: new approaches to brain analysis
Ingle et al. Advancing ASD Diagnostic Classification with Features of Continuous Wavelet Transform of fMRIand Machine Learning Algorithms
Parmar Machine learning in functional magnetic resonance neuroimaging analysis
TASSI Whole-brain analysis of resting state functional networks: a 3T fMRI twin study
Brundavanam Examining the performance of trend surface models for inference on Functional Magnetic Resonance Imaging (fMRI) data
Wang Bayesian Multimodal Local False Discovery Rate in Neuroimaging Studies
Kambli DELINEATION OF RESTING STATE FUNCTIONAL MRI INTRINSIC BRAIN NETWORKS
Gomez-Laberge et al. Selection criteria for the analysis of data-driven clusters in cerebral fMRI
Felouat Functional and Structural Analysis of fMRI/MRI Images of The Brain Using Graph Matching According to a Parallel and Distributed Approach
Lee Assessing the efficacy of pre-processing choices in fMRI: Intensity normalization effects as a function of age and task
Yu Bayesian Modeling of Complex-valued FMRI Signals
Kalus Biostatistical modeling and analysis of combined fMRI and EEG measurements
Zhao Statistical methodologies for group comparisons of brain connectivity using multimodal neuroimaging data
Welvaert Simulation of fMRI data: A statistical approach
Clement Analysis techniques for diffusion tensor imaging data
HUAIEN Functional MRI data analysis: Detection, estimation and modelling

Legal Events

Date Code Title Description
A871 Explanation of circumstances concerning accelerated examination

Free format text: JAPANESE INTERMEDIATE CODE: A871

Effective date: 20080820

A975 Report on accelerated examination

Free format text: JAPANESE INTERMEDIATE CODE: A971005

Effective date: 20080909

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20080912

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20080916

A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20080820

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20081113

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20081127

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20090113

RD03 Notification of appointment of power of attorney

Free format text: JAPANESE INTERMEDIATE CODE: A7423

Effective date: 20090210

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20090218

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20090313

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20090423

A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20090501

R150 Certificate of patent or registration of utility model

Ref document number: 4308871

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

Free format text: JAPANESE INTERMEDIATE CODE: R150

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20120515

Year of fee payment: 3

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20120515

Year of fee payment: 3

RD04 Notification of resignation of power of attorney

Free format text: JAPANESE INTERMEDIATE CODE: R3D04

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20120515

Year of fee payment: 3

S201 Request for registration of exclusive licence

Free format text: JAPANESE INTERMEDIATE CODE: R314201

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20120515

Year of fee payment: 3

R350 Written notification of registration of transfer

Free format text: JAPANESE INTERMEDIATE CODE: R350

S531 Written request for registration of change of domicile

Free format text: JAPANESE INTERMEDIATE CODE: R313531

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20120515

Year of fee payment: 3

R350 Written notification of registration of transfer

Free format text: JAPANESE INTERMEDIATE CODE: R350

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20130515

Year of fee payment: 4

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20130515

Year of fee payment: 4

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250