JP5785533B2 - 脳内電流の算出方法、算出装置およびコンピュータプログラム - Google Patents

Description

本発明は、脳内電流の算出方法、算出装置およびコンピュータプログラムに関し、特に、脳内電流分布の解像度が高いとともに、短い時間で解を得ることができる脳内電流の算出方法、算出装置およびコンピュータプログラムに関する。

近年、ロボットその他の装置を操作するのに、人の脳波（頭皮の電位）または脳磁（頭部周囲の磁界）を検出することで、人が動作を考えるだけで、ロボットその他の装置がその動作をするように制御する技術（ＢＭＩ：Ｂｒａｉｎ Ｍａｃｈｉｎｅ Ｉｎｔｅｒｆａｃｅ）が研究されている。ＢＭＩにおいては、脳波または脳磁そのものに基づいて人の脳の活動を推定することが行われることもあるが、これらの電磁的情報に基づいて、より詳細に、脳内（表面だけではなく、３次元的な脳の中）のどの部分において電流が流れたかを推定し、より正確な脳内活動を推定しようとすることも試みられている。脳波をもとに脳内電流を推定する技術はＥＥＧ逆問題と呼ばれ、脳磁をもとに推定する技術はＭＥＧ逆問題と呼ばれている。

ＥＥＧ逆問題やＭＥＧ逆問題は、電流源（脳内電流）から電位または磁場を導くリードフィールド行列を求め、このリードフィールド行列から電流源を推定することができる（例えば、特許文献１）。

ＥＥＧ逆問題やＭＥＧ逆問題を解く場合には、人体頭部を球体近似することで高速で脳内電流を計算することができるが（非特許文献１，２）、この手法は、実際の頭部形状を反映できず、正確さに問題がある。また、ＭＥＧ逆問題の場合、均質球体では半径方向を向いた電流源による外部磁界が０になるという特性があり、解析できる脳内電流が大きく制約を受けてしまうという問題がある。

そのため、実際の人体頭部形状を考慮した数値解析法を用いて電流源を計算することが試みられている（非特許文献３）。この場合、人体頭部をメッシュモデルで表現し、ＢＥＭ（Ｂｏｕｎｄａｒｙ Ｅｌｅｍｅｎｔ Ｍｅｔｈｏｄ，境界要素法）やＦＥＭ（Ｆｉｎｉｔｅ Ｅｌｅｍｅｎｔ Ｍｅｔｈｏｄ，有限要素法）を用いて電流源を求める。これらの数値解析法を用いると、人体頭部形状の影響を忠実に反映させることができるので、球体近似の制約を受けることなく脳内電流を精度よく計算することが可能になる。

特開２００３−３８４５５号公報

Sravas, "Basic mathematical and electromagnetic concepts of thebiomagnetic inverse problem," Phys. Med. Biol., vol.32, pp.11-22, 1987. D. Yao, "Electric Potential Produced by a Dipole in a HomogeneousConducting Sphere," IEEE Trans. Biomed. Eng. vol.47, pp.964-966, 2000. A. Crouzeix, "Methodes de localisation des generateurs de l’activiteelectrique cerebrale a partir de signaux electro-etmagneto-encepahlographiques," Rapport de these, Institution National desSciences Appliquees de Lyon, pp.11-31, 2001.

しかしながら、上記のような数値解析法は大量の計算が必要になるため、計算時間が長くなるという問題がある。数値解析法における計算量は、グリッドを構成する格子点の個数に大きく依存するので、グリッドのピッチを粗くすると計算量を少なくすることができるが、同時に計算精度が悪化するため、計算量を抑えつつ高精度（高分解能）な解を求めることができなかった。

本発明は、上述の背景からなされたもので、解像度が高いとともに、短い時間で脳内電流分布を得ることができる脳内電流の算出方法、算出装置およびコンピュータプログラムを提供することを目的とする。

前記した課題を解決する本発明は、コンピュータを用いて、頭部表面の電磁的情報に基づいて、脳内の電流分布を算出する脳内電流の算出方法であって、所定ピッチで配置されたグリッドを構成する格子点を設定する初期グリッド設定ステップと、グリッドを構成する格子点における電流分布を与える電流源ベクトルを［ｓ］、前記電磁的情報の配列からなる電磁ベクトルを［ｂ］、リードフィールド行列を［Ｌ］として、順問題を解くことにより［Ｌ］を求め、逆問題を解くことで［ｓ］を求め、これにより、各格子点における電流を、前記電磁的情報に基づき算出する電流算出ステップと、前記電流算出ステップで算出した各格子点における電流値に基づき、前回設定されたグリッドの一部の領域のみにおいて当該グリッドのピッチより小さなピッチでサブグリッドを構成する格子点を設定するサブグリッド設定ステップと、を備え、前記初期グリッド設定ステップと前記電流算出ステップを実行した後、前記サブグリッド設定ステップと前記電流算出ステップを１回以上繰り返すことを特徴とする。

この方法によると、初期グリッド設定ステップと電流算出ステップの実行により初期のグリッドにおける各格子点の電流値を求めることができる。そして、サブグリッド設定ステップで、求められた各格子点の電流値に基づき、すでに設定されたグリッドの一部の領域のみにおいて、そのすでに設定されたグリッドより小さなピッチでサブグリッドを構成する格子点が設定される。さらに、サブグリッドを含めたグリッドで再度電流算出ステップを行うことにより、各格子点における電流を求めることができる。

これによれば、小さなピッチで格子点が設定されるのは一部の領域に限られ、他の領域の格子点のピッチは粗いままなので、格子点の総数が少なく抑えられることで計算量を少なくし、短い時間で解を得ることができる。また、電流分布を知りたい領域においては、小さなピッチで格子点が設定されるので、必要な領域において高解像度の電流分布を得ることができる。

前記した方法は、周囲の格子点に対して電流値が高いピーク電流値を有するピーク格子点を少なくとも１つ検出するピーク格子点検出ステップをさらに備え、前記サブグリッド設定ステップにおいて、前記ピーク格子点の周囲に前記サブグリッドを構成する格子点を設定することができる。

この方法によると、電流値のピークを有する付近の領域について、電流分布を高解像度で求めることができる。

前記した方法において、前記電流算出ステップは、前回設定されたグリッド全体に対するリードフィールド行列を［Ｌ_ｉ］、今回設定されたサブグリッドの格子点（新たに設定された格子点）に対するリードフィールド行列を［Ｌ_ｓｉ］、今回設定されたグリッド全体に対するリードフィールド行列を［Ｌ_ｉ＋１］として、［Ｌ_ｉ］を［Ｌ_ｉ＋１］に対応した配列の大きさにするため、要素０の行列を追加して［Ｌ_ｉ ０］を設定し、［Ｌ_ｓｉ］を［Ｌ_ｉ＋１］に対応した配列の大きさにするため、要素０の行列を追加して［０ Ｌ_ｓｉ］を設定し、電磁ベクトル［ｂ］に基づき順問題を解くことで［０ Ｌ_ｓｉ］を求め、［Ｌ_ｉ＋１］を［Ｌ_ｉ＋１］＝［Ｌ_ｉ ０］＋［０ Ｌ_ｓｉ］により算出することが望ましい。

この方法によれば、増やした格子点に応じた要素数のリードフィールド行列［０ Ｌ_ｓｉ］のみについて順問題を解き、すでに計算済みのリードフィールド行列［Ｌ_ｉ ０］との和をとることで［Ｌ_ｉ Ｌ_ｓｉ］を得ることができるので、全体のリードフィールド行列［Ｌ_ｉ Ｌ_ｓｉ］の順問題を解く場合に比較して大幅に計算量を少なくすることができる。

前記した方法において、前記電流算出ステップは、グラム行列［Ｇ］（＝［Ｌ］［Ｌ］^Ｔ）を用いることで逆問題を解くステップを有し、前回設定されたグリッド全体に対するグラム行列を［Ｇ_ｉ］、今回設定されたグリッド全体に対するグラム行列を［Ｇ_ｉ＋１］、今回設定されたサブグリッドに対するリードフィールド行列を［Ｌ_ｓｉ］として、［Ｇ_ｉ＋１］を［Ｇ_ｉ＋１］＝［Ｇ_ｉ］＋［Ｌ_ｓｉ］［Ｌ_ｓｉ］^Ｔにより算出することが望ましい。

この方法によれば、今回設定するグリッド全体に対するグラム行列［Ｇ_ｉ＋１］を［Ｌ_ｉ＋１］［Ｌ_ｉ＋１］^Ｔを解いて求める場合に比較して計算量を少なくすることができる。

前記した課題を解決する本発明は、脳内電流の算出装置として構成することができる。すなわち、本発明は、頭部表面の電磁的情報に基づいて、脳内の電流分布を算出するための脳内電流の算出装置であって、所定ピッチで配置されたグリッドを構成する格子点を設定する初期グリッド設定手段と、グリッドを構成する格子点における電流分布を与える電流源ベクトルを［ｓ］、前記電磁的情報の配列からなる電磁ベクトルを［ｂ］、リードフィールド行列を［Ｌ］として、順問題を解くことにより［Ｌ］を求め、逆問題を解くことで［ｓ］を求め、これにより、各格子点における電流を、前記電磁的情報に基づき算出する電流算出手段と、前記電流算出手段で算出した各格子点における電流値に基づき、前回設定されたグリッドの一部の領域のみにおいて当該グリッドのピッチより小さなピッチでサブグリッドを構成する格子点を設定するサブグリッド設定手段と、前記初期グリッド設定手段と前記電流算出手段により初期のグリッドに対応した電流源ベクトル［ｓ_１］の計算を実行した後、前記サブグリッド設定手段と前記電流算出手段によるサブグリッドの設定と電流源ベクトル［ｓ_ｉ＋１］の計算を１回以上繰り返す計算実行手段と、を備える構成とすることができる。

また、前記した課題を解決する本発明は、脳内の電流分布を算出するためのコンピュータプログラムとして構成することができる。すなわち、本発明は、コンピュータを用いて、頭部表面の電磁的情報に基づいて、脳内の電流分布を算出するためのコンピュータプログラムであって、コンピュータを、所定ピッチで配置されたグリッドを構成する格子点を設定する初期グリッド設定手段、グリッドを構成する格子点における電流分布を与える電流源ベクトルを［ｓ］、前記電磁的情報の配列からなる電磁ベクトルを［ｂ］、リードフィールド行列を［Ｌ］として、順問題を解くことにより［Ｌ］を求め、逆問題を解くことで［ｓ］を求め、これにより、各格子点における電流を、前記電磁的情報に基づき算出する電流算出手段、前記電流算出手段で算出した各格子点における電流値に基づき、前回設定されたグリッドの一部の領域のみにおいて当該グリッドのピッチより小さなピッチでサブグリッドを構成する格子点を設定するサブグリッド設定手段、前記初期グリッド設定手段と前記電流算出手段により初期のグリッドに対応した電流源ベクトル［ｓ_１］の計算を実行した後、前記サブグリッド設定手段と前記電流算出手段によるサブグリッドの設定と電流源ベクトル［ｓ_ｉ＋１］の計算を１回以上繰り返す計算実行手段、として機能させる構成とすることができる。

本発明の脳内電流の算出方法、算出装置およびコンピュータプログラムによれば、高い解像度で脳内電流分布を求めることができるとともに、短い時間で解を得ることができる。

実施形態に係る脳内電流の算出装置のブロック図である。 サブグリッドの設定を説明する図であり、（ａ）設定前と（ｂ）設定後である。 脳内電流の算出処理を説明するフローチャートである。 頭部の表面メッシュモデルと脳内の初期グリッドの一例を示す図である。 脳内に１回目のサブグリッドを設定したときの頭部の表面メッシュモデルと脳内グリッドの一例を示す図である。 脳内に２回目のサブグリッドを設定したときの頭部の表面メッシュモデルと脳内グリッドの一例を示す図である。

次に、本発明の一実施形態について、適宜図面を参照しながら詳細に説明する。一実施形態に係る脳内電流の算出装置１００は、頭皮の表面の電位を測定することで得られる脳波や、頭皮の近傍の磁界を検出することで得られる脳磁などの電磁ベクトル［ｂ］（脳波・脳磁の情報の集合をベクトル表現したもの）に基づき、脳内の電流分布、具体的には、電磁ベクトル［ｂ］を発生した原因である電流源ベクトル［ｓ］を算出する装置である。なお、電磁ベクトル［ｂ］の要素の数は、センサの数と同じであるので、通常、十数個〜数十個であり、電流源ベクトル［ｓ］の要素の数は、脳内モデルのグリッドの格子点の数と同じであるので数千個〜数万個である。

上記の計算のため、算出装置１００は、図１に示すように、脳内電流を算出するための手段として、初期グリッド設定手段１１０、電流算出手段１２０、ピーク格子点検出手段１３０、サブグリッド設定手段１４０および計算実行手段１５０を有している。算出装置１００は、図示しないＣＰＵ、ＲＯＭ、ＲＡＭおよび外部記憶装置などの記憶装置を有するコンピュータからなり、ＲＡＭにコンピュータプログラムが適宜読み込まれてＣＰＵで実行されることでこれらの各手段が実現されている。また、算出装置１００は、グリッドの座標データや、リードフィールド行列［Ｌ］、電磁ベクトル［ｂ］、電流源ベクトル［ｓ］、グラム行列［Ｇ］などの各種の変数などを記憶するための記憶部１９０を有している。

初期グリッド設定手段１１０は、所定ピッチで配置されたグリッドを構成する格子点である初期グリッドを設定する手段である。この初期グリッドの設定は、公知の手法と同様に空間座標を定義することで行われる。初期グリッドのピッチは、本発明を実施しようとするコンピュータの計算能力に応じて、比較的短時間で電流の算出が行えるような大きな（粗い）ものとするとよい。例えば、脳の実寸法に対応させると１０ｍｍ程度のピッチで設定することができ、この場合、１２００個程度の格子点数に収めることができる。

電流算出手段１２０は、初期グリッド設定手段１１０およびサブグリッド設定手段１４０によりすでに設定されたグリッドを用いて、電磁ベクトル［ｂ］から電流源ベクトル［ｓ］を算出する手段である。ここでのグリッドは、初期グリッドとサブグリッドの両方を含む。電流を算出するための基本的な考え方は、グリッドを構成する格子点における電流分布を与える電流源ベクトルを［ｓ］、電磁的情報の配列からなる電磁ベクトルを［ｂ］、リードフィールド行列を［Ｌ］として、順問題を解くことにより［Ｌ］を求め、逆問題を解くことで［ｓ］を求め、これにより、各格子点における電流を、前記電磁的情報に基づき算出するというものであり、この点においては従来と同様である。

少し詳しく説明すれば、電流源ベクトル［ｓ］、リードフィールド行列［Ｌ］および電磁ベクトル［ｂ］は、
［Ｌ］［ｓ］＝［ｂ］
の関係にあり、まず、順問題を解く際には、［ｓ］のうちｋ番目の電流成分を１とし、他の電流成分を０と置いて［ｂ］を求める。得られた［ｂ］ベクトルがリードフィールド行列［Ｌ］のｋ番目の列ベクトル成分となる。この手順を順に繰り返して［Ｌ］の要素を決定していく。
そして、リードフィールド行列［Ｌ］が決定した後は、上記式の逆問題を解く各種の手法により電流源ベクトル［ｓ］を解くことができる。

逆問題を解くには、例えば、以下の各手法がある。なお、電流源ベクトル［ｓ］は、前記したように電磁ベクトル［ｂ］より要素数が多く未知数に対して方程式の数が少ないため、一意には解が定まらない。以下の電流源ベクトル［ｓ］はすべて推定解であり、便宜上、同じ［ｓ］により表記する。
電流源ベクトル［ｓ］を解く式を重み行列［Ｗ］を用いて一般式で
［ｓ］＝［Ｗ］^Ｔ［ｂ］
で表すと、［Ｗ］は、例えば以下（１）から（１０）の各手法で解くことができる。なお、以下の各式において［Ｇ］はグラム行列（＝［Ｌ］［Ｌ］^Ｔ）であり、［Ｒ］は［ｂ］の共分散行列である。

（１）Ｍｉｎｉｍｕｍ Ｎｏｒｍ法
［Ｗ］＝［Ｇ］^−１［Ｌ］

（２）重み付Ｍｉｎｉｍｕｍ Ｎｏｒｍ法
［Ｗ］＝［Ｇ］^−１［Ｌ］（［Ｌ］^Ｔ［Ｇ］^−２［Ｌ］）^−１／２

（３）ＷＲＯＰ法(Weighted Resolution Optimization)
［Ｗ］＝［Ｇ］^−１［Ｌ］（［Ｌ］^Ｔ［Ｇ］^−１［Ｌ］）^−１

（４）重み付ＷＲＯＰ法

（５）ｓＬＯＲＥＴＡ法(standardized LOw REsolution brain electromagnetic Tomography)
［Ｗ］＝［Ｇ］^−１［Ｌ］（［Ｌ］^Ｔ［Ｇ］^−１［Ｌ］）^−１／２

（６）Ｍｉｎｉｍｕｍ Ｖａｒｉａｎｃｅ版Ｍｉｎｉｍｕｍ Ｎｏｒｍ法
［Ｗ］＝［Ｒ］^−１［Ｌ］

（７）重み付Ｍｉｎｉｍｕｍ Ｖａｒｉａｎｃｅ法
［Ｗ］＝［Ｒ］^−１［Ｌ］（［Ｌ］^Ｔ［Ｒ］^−２［Ｌ］）^−１／２

（８）Ｍｉｎｉｍｕｍ Ｖａｒｉａｎｃｅ法
［Ｗ］＝［Ｒ］^−１［Ｌ］（［Ｌ］^Ｔ［Ｒ］^−１［Ｌ］）^−１

（９）Ａｒｒａｙ−ｇａｉｎ ｃｏｎｓｔｒａｉｎｔ Ｍｉｎｉｍｕｍ Ｖａｒｉａｎｃｅ法

（１０）Ｍｉｎｉｍｕｍ Ｖａｒｉａｎｃｅ版ｓＬＯＲＥＴＡ法
［Ｗ］＝［Ｒ］^−１［Ｌ］（［Ｌ］^Ｔ［Ｒ］^−１［Ｌ］）^−１／２

そして、電流算出手段１２０は、サブグリッド設定手段１４０によりサブグリッドが設定された、２回目以降の計算においては、今回の（つまり、最後にサブグリッドが設定された後の）グリッド全体に対するリードフィールド行列［Ｌ_ｉ＋１］を次のようにして求める。
まず、前回設定されたグリッド全体に対するリードフィールド行列を［Ｌ_ｉ］、今回設定されたサブグリッドに対するリードフィールド行列を［Ｌ_ｓｉ］として、［Ｌ_ｉ］を［Ｌ_ｉ＋１］に対応した配列の大きさにする。ここで、［Ｌ_ｉ＋１］は、［Ｌ_ｉ］に対し、行数は同じであるが、列数が追加したサブグリッドの未知数に対応して増えているので、［Ｌ_ｉ］に要素０の行列を追加して［Ｌ_ｉ ０］を設定する。同様に、［Ｌ_ｓｉ］を［Ｌ_ｉ＋１］に対応した配列の大きさにするため、要素０の行列を追加して［０ Ｌ_ｓｉ］を設定する。そして、入力された電磁ベクトル［ｂ］に基づき順問題を解くことで［０ Ｌ_ｓｉ］を求め、［Ｌ_ｉ＋１］を
［Ｌ_ｉ＋１］＝［Ｌ_ｉ ０］＋［０ Ｌ_ｓｉ］
により算出する。

これにより、サブグリッド設定手段１４０により増やした格子点に応じた要素数のリードフィールド行列［０ Ｌ_ｓｉ］のみについて順問題を解き、すでに計算済みのリードフィールド行列［Ｌ_ｉ ０］との和をとることで［Ｌ_ｉ Ｌ_ｓｉ］（＝［Ｌ_ｉ＋１］）を得ることができる。電流源ベクトル［ｓ］を解くまでの計算の過程で、順問題を解くための計算の量は大きな割合を占めるので、全体のリードフィールド行列［Ｌ_ｉ Ｌ_ｓｉ］の順問題を解く場合に比較して大幅に計算量を少なくすることができる。

また、電流算出手段１２０は、前記した逆問題の解法において、グラム行列［Ｇ］を用いる場合には、前回設定されたグリッド全体に対するグラム行列を［Ｇ_ｉ］、今回設定されたグリッド全体に対するグラム行列を［Ｇ_ｉ＋１］、今回設定されたサブグリッドに対するリードフィールド行列を［Ｌ_ｓｉ］として、［Ｇ_ｉ＋１］を
［Ｇ_ｉ＋１］＝［Ｇ_ｉ］＋［Ｌ_ｓｉ］［Ｌ_ｓｉ］^Ｔ
により算出することが望ましい。
このように構成しても、以下のように［Ｌ_ｉ＋１］［Ｌ_ｉ＋１］^Ｔを計算する場合と同じ解を得ることができる。

［Ｇ_ｉ＋１］＝［Ｌ_ｉ＋１］［Ｌ_ｉ＋１］^Ｔ
＝［［Ｌ_ｉ ０］＋［０ Ｌ_ｓｉ］］［［Ｌ_ｉ ０］＋［０ Ｌ_ｓｉ］］^Ｔ
＝［Ｌ_ｉ ０］［Ｌ_ｉ ０］^Ｔ＋［Ｌ_ｉ ０］［０ Ｌ_ｓｉ］^Ｔ
＋［０ Ｌ_ｓｉ］［Ｌ_ｉ ０］^Ｔ＋［０ Ｌ_ｓｉ］［０ Ｌ_ｓｉ］^Ｔ
＝［Ｌ_ｉ］［Ｌ_ｉ］^Ｔ＋［０］＋［０］＋［Ｌ_ｓｉ］［Ｌ_ｓｉ］^Ｔ
＝［Ｇ_ｉ］＋［Ｌ_ｓｉ］［Ｌ_ｓｉ］^Ｔ

このように電流算出手段１２０を構成すると、［Ｇ_ｉ］はすでに計算済みであるので、新たに［Ｌ_ｓｉ］［Ｌ_ｓｉ］^Ｔを計算して［Ｇ_ｉ］との和をとることで、少ない計算量で［Ｇ_ｉ＋１］を得ることができる。

ピーク格子点検出手段１３０は、電流算出手段１２０により算出された、各格子点の電流に基づき、周囲の格子点に対して電流値が高いピーク電流値を有するピーク格子点を少なくとも１つ検出する手段である。ピーク検出の方法は、公知の各種の方法を用いることができるが、例えば、上下、左右、前後に隣接するいずれの格子点よりも高い電流値を有する格子点を複数検出し、検出された格子点の中から上位の所定個数の格子点を選択することでピーク電流値を有する格子点を検出することができる。検出するピーク格子点の個数は、計算対象となるデータ（脳波、脳磁）に応じて変えることができる。あまり多くのピーク格子点を選択すると、細かいノイズまで選択してしまうことで、サブグリッドを広い範囲に設定することになり、計算量が多くなるので、ピーク格子点の数は、逆問題を解く手法や予想される脳内電流のピーク数に応じて適宜決定されるが、通常は５〜５０個程度の範囲が良く、できるだけ少ない方が望ましい。ピークの検出が不十分な場合は、ピーク格子点の数を増やすと良い。

サブグリッド設定手段１４０は、電流算出手段１２０が算出した各格子点における電流値に基づき、前回設定されたグリッドの一部の領域のみにおいて当該グリッドのピッチより小さなピッチでサブグリッドを構成する格子点を設定する手段である。本実施形態においては、ピーク格子点付近の電流分布を高解像度で計算するため、ピーク格子点検出手段１３０が検出したピーク格子点付近にサブグリッドを設定する。

例えば、サブグリッド設定手段１４０は、図２（ａ）に示すように、ピーク格子点ＧＰの周囲にすでにピッチをＰ_ｉの格子点Ｇ_ｉからなるグリッドＧが設定されている場合に、図２（ｂ）に示すように、ピーク格子点ＧＰを中心とする、上下左右前後にＮｒ×Ｐ_ｉ以内の立方体（図２は、便宜上２次元で表現している。）の範囲でピッチＰ_ｉより小さなピッチＰ_ｉ＋１のグリッドができるように新たな格子点Ｇ_ｉ＋１を設定する。本発明では、この格子点Ｇ_ｉを補完する新たな格子点Ｇ_ｉ＋１を設定することでできたピッチＰ_ｉ＋１のグリッドの部分をサブグリッドＳＧ_ｉ＋１と呼ぶ。図２（ｂ）では、Ｎｒ＝３とし、ピッチＰ_ｉ＋１をＰ_ｉ／２としているが、これらは任意であり、例えば、Ｎｒは２〜５程度とすることができるし、ピッチＰ_ｉ＋１は、Ｐ_ｉ／２〜Ｐ_ｉ／５程度とすることができる。また、サブグリッドＳＧ_ｉ＋１を設定する範囲も、立方体に限らず、球体とすることができる。さらに、ピーク格子点ＧＰが線状に連なっている場合には、サブグリッドＳＧ_ｉ＋１を設定する範囲をその線に沿って延びる直方体または円柱とし、線から直交する方向のＮｒを小さくすることで、サブグリッドＳＧ_ｉ＋１の大きさを小さくしつつ、ピーク格子点ＧＰ付近の電流分布の解像度を高くすることができる。

計算実行手段１５０は、初期グリッド設定手段１１０による初期グリッドの設定と、初期グリッドに対する電流算出手段１２０による電流源ベクトル［ｓ_１］の算出を実行した後、サブグリッド設定手段１４０によるサブグリッドの設定と、サブグリッドが設定された新たなグリッドについて電流算出手段１２０による電流源ベクトル［ｓ_ｉ＋１］の算出を１回以上繰り返す手段である。すなわち、計算実行手段１５０は、コンピュータプログラムにおけるメインフローに相当する処理手段である。

計算実行手段１５０によるサブグリッド設定ステップと電流算出ステップの繰り返し数は任意であり、計算対象とするデータに応じて設定するとよい。

以上のような構成の算出装置１００による脳内電流の算出処理（算出方法）を図３のフローチャートおよび図４〜図６のモデルを参照しながら説明する。
図３の処理の全体は、計算実行手段１５０により制御され、まず、初期グリッド設定手段１１０により脳内に相当する部分に、ピッチＰ_１の初期グリッドＩＧが設定される（Ｓ１、初期グリッド設定ステップ、図４参照）。

そして、電流算出手段１２０は、初期グリッドＩＧについて、入力された電磁ベクトル［ｂ］に基づいて順問題を解くことによりリードフィールド行列［Ｌ_１］を算出する（Ｓ２、電流算出ステップ）。そして、逆問題を解く手法に応じて、全領域についての他の変数、例えば、グラム行列［Ｇ］や共分散行列［Ｒ］、重み行列［Ｗ］などを算出する。算出した各変数は、記憶部１９０に記憶される。なお、以下においては、記憶部１９０への算出結果や算出過程の記憶についての説明は省略する。

そして、重み行列［Ｗ］を用いて逆問題を解くことにより、各格子点Ｇ_１の電流値（つまり、電流源ベクトル［ｓ_１］）を算出する（Ｓ３、電流算出ステップ）。電流値の算出の後、繰り返し計算実行手段１５０は、反復回数の上限値か否か判定し、上限値に達していれば（Ｓ４，Ｙｅｓ）、処理を終了し、上限値に達していなければ（Ｓ４，Ｎｏ）、サブグリッドの設定と電流値の再計算のためにステップＳ５に進む。

反復回数の上限値に達していない場合、ピーク格子点検出手段１３０は、ピーク値を検出する（Ｓ５、ピーク格子点検出ステップ）。これにより、例えば、図４の例においては、表示された断面において３つのピーク格子点ＧＰが検出されている。

次に、サブグリッド設定手段１４０は、ピーク格子点検出手段１３０が検出したピーク格子点ＧＰに基づいて、ピーク格子点ＧＰの周囲の所定の範囲に新たな格子点Ｇ_２を設定することでサブグリッドＳＧ_２を設定する（Ｓ６、サブグリッド設定ステップ、図５参照）。そして、サブグリッド領域について、変数を計算する（Ｓ７、電流算出ステップ）。すなわち、新たに設定した格子点Ｇ_２に対応したリードフィールド行列［Ｌ_ｓ１］を［Ｌ_２］に対応した配列の大きさにするため、要素０の行列を追加して［０ Ｌ_ｓ１］を設定する。そして、電磁ベクトル［ｂ］に基づき順問題を解くことで［０ Ｌ_ｓ１］を求める。

次に、グリッド全体の変数を計算して更新する（Ｓ８、電流算出ステップ）。ここでの変数計算は、リードフィールド行列［Ｌ_２］について、前記したように、［Ｌ_１］を［Ｌ_２］に対応した配列の大きさにするため、要素０の行列を追加して［Ｌ_１ ０］を設定し、新たに設定した［Ｌ_２］を
［Ｌ_２］＝［Ｌ_１ ０］＋［０ Ｌ_ｓ１］
により算出する。また、逆問題を解くのにグラム行列［Ｇ］を使用する場合には、
［Ｇ_２］＝［Ｇ_１］＋［Ｌ_ｓ１］［Ｌ_ｓ１］^Ｔ
により算出する。

その後、計算実行手段１５０は、ステップＳ３に戻って電流値の計算を行い、必要な回数、サブグリッドＳＧ_ｉの設定と電流源ベクトル［ｓ_ｉ］の計算を繰り返す（Ｓ４〜Ｓ８，Ｓ３）。
このようにして、例えば、図４〜図６の例の場合、図５に示すようにサブグリッドＳＧ_２の設定後、再度ピーク格子点ＧＰを検出し、図６に示すようにその周囲にさらに細かいサブグリッドＳＧ_３を設定する。図６に示すように、サブグリッドＳＧ_３の設定後に電流源ベクトル［ｓ_３］を計算することで、ピーク格子点ＧＰの周辺において高解像度な電流分布を得ることができる。

以上のように、本実施形態の脳内電流の算出装置１００（コンピュータプログラム）、および算出装置１００による算出方法によれば、ピーク格子点の周辺には、細かいサブグリッドを設定することで、高解像度の電流分布を得ることができるとともに、ピーク格子点から遠い部分は、粗いピッチのグリッドが設定されることで、格子点の総数が大きくなりすぎないので、計算量を少なくして短時間で脳内電流分布を算出することができる。

また、２回目以降の電流値の算出において、リードフィールド行列［Ｌ_ｉ＋１］を
［Ｌ_ｉ＋１］＝［Ｌ_ｉ ０］＋［０ Ｌ_ｓｉ］
により算出することで、計算量を少なくして短時間で脳内電流分布を得ることができる。また、逆問題を解くのにグラム行列［Ｇ］を使用する場合には、グラム行列［Ｇ_ｉ＋１］を
［Ｇ_ｉ＋１］＝［Ｇ_ｉ］＋［Ｌ_ｓｉ］［Ｌ_ｓｉ］^Ｔ
により算出することで計算量を少なくして短時間で脳内電流分布を得ることができる。

この計算量の比較について、簡単に説明する。
基準のグリッド（例えば、初期グリッド）のピッチで脳内電流分布を計算するのに掛かる時間をｔとすると、全体を１／２のピッチのグリッドにした場合には８ｔ、全体を１／４のピッチのグリッドにした場合には６４ｔの時間がかかる。

これに対し、本実施形態のように、一部にのみサブグリッドを形成して計算する場合、ピーク格子点の周囲にサブグリッドを設定する範囲Ｎｒ＝２、ピーク格子点の個数Ｎｐ＝５の条件下で、繰り返し計算回数を１回とし、サブグリッドのピッチを初期グリッドの１／２にした場合には、計算時間は２ｔ＋ｇとなる。また、繰り返し計算回数を２回とし、２回目のサブグリッドのピッチを初期グリッドの１／４にした場合には、計算時間は４ｔ＋２ｇとなる。

なお、ｇは、下記の（１）〜（３）のような演算処理に要する時間である。
（１）サブグリッドを生成することでリードフィールド行列のサイズが大きくなるため、元のリードフィールド行列に対して増加したデータが格納できるように計算領域を拡張（メモリー確保）し、この領域にサブグリッドのデータを書き込むための時間。
（２）サブグリッドを追加したことにより、計算結果を表示するときの六面体要素が変更されるが、この新たな六面体要素を構成する格子点の情報を書き換える時間。
（３）検出されたピーク格子点が互いに近い場合、これらの周囲のサブグリッドを構成する格子点が重複する場合があるので、この重複した格子点を１つにまとめるための処理に必要な時間。
これらの処理に必要な時間であるｇは、ｔに比べて遙かに短い時間であるので、本実施形態の算出装置１００によれば非常に短い時間で脳内電流分布を得ることができることが分かる。

以上に本発明の実施形態について説明したが、本発明は、前記した実施形態に限定されることなく、適宜変形して実施することができる。
前記実施形態においては、サブグリッドを設定する領域を、ピーク格子点の周囲としたが、サブグリッドを設定する領域は、ピーク格子点の周囲には限定されず、例えば、所定の電流値を有する格子点の周囲とすることもできる。

１００ 算出装置
１１０ 初期グリッド設定手段
１２０ 電流算出手段
１３０ ピーク格子点検出手段
１４０ サブグリッド設定手段
１５０ 計算実行手段
１９０ 記憶部

Claims (7)

1. コンピュータを用いて、頭部表面の電磁的情報に基づいて、脳内の電流分布を算出する脳内電流の算出方法であって、
   所定ピッチで配置されたグリッドを構成する格子点を設定する初期グリッド設定ステップと、
   グリッドを構成する格子点における電流分布を与える電流源ベクトルを［ｓ］、前記電磁的情報の配列からなる電磁ベクトルを［ｂ］、リードフィールド行列を［Ｌ］として、順問題を解くことにより［Ｌ］を求め、逆問題を解くことで［ｓ］を求め、これにより、各格子点における電流を、前記電磁的情報に基づき算出する電流算出ステップと、
   前記電流算出ステップで算出した各格子点における電流値に基づき、前回設定されたグリッドの一部の領域のみにおいて当該グリッドのピッチより小さなピッチでサブグリッドを構成する格子点を設定するサブグリッド設定ステップと、を備え、
   前記電流算出ステップは、
   前回設定されたグリッド全体に対するリードフィールド行列を［Ｌ _ｉ ］、今回設定されたサブグリッドに対するリードフィールド行列を［Ｌ _ｓｉ ］、今回設定されたグリッド全体に対するリードフィールド行列を［Ｌ _ｉ＋１ ］として、［Ｌ _ｉ ］を［Ｌ _ｉ＋１ ］に対応した配列の大きさにするため、要素０の行列を追加して［Ｌ _ｉ ０］を設定し、［Ｌ _ｓｉ ］を［Ｌ _ｉ＋１ ］に対応した配列の大きさにするため、要素０の行列を追加して［０ Ｌ _ｓｉ ］を設定し、
   電磁ベクトル［ｂ］に基づき順問題を解くことで［０ Ｌ _ｓｉ ］を求め、［Ｌ _ｉ＋１ ］を
   ［Ｌ _ｉ＋１ ］＝［Ｌ _ｉ ０］＋［０ Ｌ _ｓｉ ］
   により算出し、
   前記初期グリッド設定ステップと前記電流算出ステップを実行した後、前記サブグリッド設定ステップと前記電流算出ステップを１回以上繰り返すことを特徴とする脳内電流の算出方法。
2. コンピュータを用いて、頭部表面の電磁的情報に基づいて、脳内の電流分布を算出する脳内電流の算出方法であって、
   所定ピッチで配置されたグリッドを構成する格子点を設定する初期グリッド設定ステップと、
   グリッドを構成する格子点における電流分布を与える電流源ベクトルを［ｓ］、前記電磁的情報の配列からなる電磁ベクトルを［ｂ］、リードフィールド行列を［Ｌ］として、順問題を解くことにより［Ｌ］を求め、逆問題を解くことで［ｓ］を求め、これにより、各格子点における電流を、前記電磁的情報に基づき算出する電流算出ステップと、
   前記電流算出ステップで算出した各格子点における電流値に基づき、前回設定されたグリッドの一部の領域のみにおいて当該グリッドのピッチより小さなピッチでサブグリッドを構成する格子点を設定するサブグリッド設定ステップと、を備え、
   前記電流算出ステップは、
   グラム行列［Ｇ］（＝［Ｌ］［Ｌ］ ^Ｔ ）を用いることで逆問題を解くステップを有し、
   前回設定されたグリッド全体に対するグラム行列を［Ｇ _ｉ ］、今回設定されたグリッド全体に対するグラム行列を［Ｇ _ｉ＋１ ］、今回設定されたサブグリッドに対するリードフィールド行列を［Ｌ _ｓｉ ］として、［Ｇ _ｉ＋１ ］を
   ［Ｇ _ｉ＋１ ］＝［Ｇ _ｉ ］＋［Ｌ _ｓｉ ］［Ｌ _ｓｉ ］ ^Ｔ
   により算出し、
   前記初期グリッド設定ステップと前記電流算出ステップを実行した後、前記サブグリッド設定ステップと前記電流算出ステップを１回以上繰り返すことを特徴とする脳内電流の算出方法。
3. 周囲の格子点に対して電流値が高いピーク電流値を有するピーク格子点を少なくとも１つ検出するピーク格子点検出ステップをさらに備え、
   前記サブグリッド設定ステップにおいて、前記ピーク格子点の周囲に前記サブグリッドを構成する格子点を設定することを特徴とする請求項１または請求項２に記載の脳内電流の算出方法。
4. 頭部表面の電磁的情報に基づいて、脳内の電流分布を算出するための脳内電流の算出装置であって、
   所定ピッチで配置されたグリッドを構成する格子点を設定する初期グリッド設定手段と、
   グリッドを構成する格子点における電流分布を与える電流源ベクトルを［ｓ］、前記電磁的情報の配列からなる電磁ベクトルを［ｂ］、リードフィールド行列を［Ｌ］として、順問題を解くことにより［Ｌ］を求め、逆問題を解くことで［ｓ］を求め、これにより、各格子点における電流を、前記電磁的情報に基づき算出する電流算出手段と、
   前記電流算出手段で算出した各格子点における電流値に基づき、前回設定されたグリッドの一部の領域のみにおいて当該グリッドのピッチより小さなピッチでサブグリッドを構成する格子点を設定するサブグリッド設定手段と、
   前記初期グリッド設定手段と前記電流算出手段により初期のグリッドに対応した電流源ベクトル［ｓ_１］の計算を実行した後、前記サブグリッド設定手段と前記電流算出手段によるサブグリッドの設定と電流源ベクトル［ｓ_ｉ＋１］の計算を１回以上繰り返す計算実行手段と、を備え、
   前記電流算出手段は、
   前回設定されたグリッド全体に対するリードフィールド行列を［Ｌ _ｉ ］、今回設定されたサブグリッドに対するリードフィールド行列を［Ｌ _ｓｉ ］、今回設定されたグリッド全体に対するリードフィールド行列を［Ｌ _ｉ＋１ ］として、［Ｌ _ｉ ］を［Ｌ _ｉ＋１ ］に対応した配列の大きさにするため、要素０の行列を追加して［Ｌ _ｉ ０］を設定し、［Ｌ _ｓｉ ］を［Ｌ _ｉ＋１ ］に対応した配列の大きさにするため、要素０の行列を追加して［０ Ｌ _ｓｉ ］を設定し、
   電磁ベクトル［ｂ］に基づき順問題を解くことで［０ Ｌ _ｓｉ ］を求め、［Ｌ _ｉ＋１ ］を
   ［Ｌ _ｉ＋１ ］＝［Ｌ _ｉ ０］＋［０ Ｌ _ｓｉ ］
   により算出することを特徴とする脳内電流の算出装置。
5. 頭部表面の電磁的情報に基づいて、脳内の電流分布を算出するための脳内電流の算出装置であって、
   所定ピッチで配置されたグリッドを構成する格子点を設定する初期グリッド設定手段と、
   グリッドを構成する格子点における電流分布を与える電流源ベクトルを［ｓ］、前記電磁的情報の配列からなる電磁ベクトルを［ｂ］、リードフィールド行列を［Ｌ］として、順問題を解くことにより［Ｌ］を求め、逆問題を解くことで［ｓ］を求め、これにより、各格子点における電流を、前記電磁的情報に基づき算出する電流算出手段と、
   前記電流算出手段で算出した各格子点における電流値に基づき、前回設定されたグリッドの一部の領域のみにおいて当該グリッドのピッチより小さなピッチでサブグリッドを構成する格子点を設定するサブグリッド設定手段と、
   前記初期グリッド設定手段と前記電流算出手段により初期のグリッドに対応した電流源ベクトル［ｓ_１］の計算を実行した後、前記サブグリッド設定手段と前記電流算出手段によるサブグリッドの設定と電流源ベクトル［ｓ_ｉ＋１］の計算を１回以上繰り返す計算実行手段と、を備え、
   前記電流算出手段は、
   グラム行列［Ｇ］（＝［Ｌ］［Ｌ］ ^Ｔ ）を用いることで逆問題を解き、
   前回設定されたグリッド全体に対するグラム行列を［Ｇ _ｉ ］、今回設定されたグリッド全体に対するグラム行列を［Ｇ _ｉ＋１ ］、今回設定されたサブグリッドに対するリードフィールド行列を［Ｌ _ｓｉ ］として、［Ｇ _ｉ＋１ ］を
   ［Ｇ _ｉ＋１ ］＝［Ｇ _ｉ ］＋［Ｌ _ｓｉ ］［Ｌ _ｓｉ ］ ^Ｔ
   により算出することを特徴とする脳内電流の算出装置。
6. コンピュータを用いて、頭部表面の電磁的情報に基づいて、脳内の電流分布を算出するためのコンピュータプログラムであって、コンピュータを、
   所定ピッチで配置されたグリッドを構成する格子点を設定する初期グリッド設定手段、
   グリッドを構成する格子点における電流分布を与える電流源ベクトルを［ｓ］、前記電磁的情報の配列からなる電磁ベクトルを［ｂ］、リードフィールド行列を［Ｌ］として、順問題を解くことにより［Ｌ］を求め、逆問題を解くことで［ｓ］を求め、これにより、各格子点における電流を、前記電磁的情報に基づき算出する電流算出手段、
   前記電流算出手段で算出した各格子点における電流値に基づき、前回設定されたグリッドの一部の領域のみにおいて当該グリッドのピッチより小さなピッチでサブグリッドを構成する格子点を設定するサブグリッド設定手段、
   前記初期グリッド設定手段と前記電流算出手段により初期のグリッドに対応した電流源ベクトル［ｓ_１］の計算を実行した後、前記サブグリッド設定手段と前記電流算出手段によるサブグリッドの設定と電流源ベクトル［ｓ_ｉ＋１］の計算を１回以上繰り返す計算実行手段、
   として機能させ、
   前記電流算出手段は、
   前回設定されたグリッド全体に対するリードフィールド行列を［Ｌ _ｉ ］、今回設定されたサブグリッドに対するリードフィールド行列を［Ｌ _ｓｉ ］、今回設定されたグリッド全体に対するリードフィールド行列を［Ｌ _ｉ＋１ ］として、［Ｌ _ｉ ］を［Ｌ _ｉ＋１ ］に対応した配列の大きさにするため、要素０の行列を追加して［Ｌ _ｉ ０］を設定し、［Ｌ _ｓｉ ］を［Ｌ _ｉ＋１ ］に対応した配列の大きさにするため、要素０の行列を追加して［０ Ｌ _ｓｉ ］を設定し、
   電磁ベクトル［ｂ］に基づき順問題を解くことで［０ Ｌ _ｓｉ ］を求め、［Ｌ _ｉ＋１ ］を
   ［Ｌ _ｉ＋１ ］＝［Ｌ _ｉ ０］＋［０ Ｌ _ｓｉ ］
   により算出することを特徴とするコンピュータプログラム。
7. コンピュータを用いて、頭部表面の電磁的情報に基づいて、脳内の電流分布を算出するためのコンピュータプログラムであって、コンピュータを、
   所定ピッチで配置されたグリッドを構成する格子点を設定する初期グリッド設定手段、
   グリッドを構成する格子点における電流分布を与える電流源ベクトルを［ｓ］、前記電磁的情報の配列からなる電磁ベクトルを［ｂ］、リードフィールド行列を［Ｌ］として、順問題を解くことにより［Ｌ］を求め、逆問題を解くことで［ｓ］を求め、これにより、各格子点における電流を、前記電磁的情報に基づき算出する電流算出手段、
   前記電流算出手段で算出した各格子点における電流値に基づき、前回設定されたグリッドの一部の領域のみにおいて当該グリッドのピッチより小さなピッチでサブグリッドを構成する格子点を設定するサブグリッド設定手段、
   前記初期グリッド設定手段と前記電流算出手段により初期のグリッドに対応した電流源ベクトル［ｓ_１］の計算を実行した後、前記サブグリッド設定手段と前記電流算出手段によるサブグリッドの設定と電流源ベクトル［ｓ_ｉ＋１］の計算を１回以上繰り返す計算実行手段、
   として機能させ、
   前記電流算出手段は、
   グラム行列［Ｇ］（＝［Ｌ］［Ｌ］ ^Ｔ ）を用いることで逆問題を解き、
   前回設定されたグリッド全体に対するグラム行列を［Ｇ _ｉ ］、今回設定されたグリッド全体に対するグラム行列を［Ｇ _ｉ＋１ ］、今回設定されたサブグリッドに対するリードフィールド行列を［Ｌ _ｓｉ ］として、［Ｇ _ｉ＋１ ］を
   ［Ｇ _ｉ＋１ ］＝［Ｇ _ｉ ］＋［Ｌ _ｓｉ ］［Ｌ _ｓｉ ］ ^Ｔ
   により算出することを特徴とするコンピュータプログラム。

Images