JP6613161B2 - Ｍｒｉシミュレーションシステム - Google Patents

Description

本発明はＭＲＩシミュレーションシステムに関する。

近年、静磁場強度の増大、勾配磁場系の高速化、ＲＦコイルの多チャンネル化、そして送受信系のデジタル化などのＭＲＩハードウェアの進化に伴い、ＭＲＩのパルスシーケンスや画像再構成手法も高度かつ複雑になってきている。このような状況から、ＭＲＩハードウェアを必要とするシーケンス開発には、多大な時間とコストを要するようになってきている。このため、ＩＤＥＡ（ＳＩＥＭＥＮＳ）やＥＰＩＣ（ＧＥ）などの統合的なＭＲＩシーケンス開発環境が提供されており、これらはシミュレータを適宜利用していると思われる。ところが、このような開発環境は、特定のメーカの機種に特化されており、機種によらない汎用的なものは提案されてはいるものの、実装には装置側の情報公開が不可欠なこともあり、実装のハードルは極めて高いものとなっている。実機ＭＲＩ装置とＭＲＩシミュレーションの相対関係の詳細は明らかにはされていない（非特許文献１、非特許文献２参照）。

ＭＲＩ数値シミュレーションの方法として、入力されたパルスシーケンス情報、装置特性情報、検査対象情報を内部データに変換し、検査対象情報生成部は、検査対象情報の磁化分布を格子に分割し、各格子に分布に応じた大きさ、緩和時間を持つ磁化を配置し、装置特性情報生成部は、検査対象情報生成部と同じ大きさの格子に視野を分割し、各格子での磁場強度を求めて内部データとし、磁化グループ生成部は、検査対象を等磁場領域に分割し、磁化の置かれている位置の磁場強度が等しいものを同グループとして、磁化を複数のグループに分け、結果をメモリに格納する方法が知られている（特許文献１参照）。

Magland JF, Li C, Langham MC, WehrliFW,Pulse sequence programming in a dynamic visual environment: SequenceTree. MagnReson Med Mar 7(2015) K.Layton, S.Kroboth, J.Leupold, H.Yu, F.Jia, S.Littin, T.Stocker, M.Zaitsev. A hardware-independent environment for MR acquisition and simulation. 23rd ISMRM annualmeeting(Toronto), 2015, p2484

特開平０９−０４７４４２号公報

ＭＲＩシミュレーションは、ＮＭＲ現象を取り扱うので、単純なＮＭＲシミュレーションを包含している。工夫をすれば汎用の計算機を用いて、ＮＭＲ現象やＭＲＩ現象のシミュレーションが実現できる状況となっている。一方、ＮＭＲ／ＭＲＩ現象を実験的に取り扱う場合、緩和時間という避けられない待ち時間があるため、ＭＲＩ装置を複数用意するということしか、開発速度を上げることができなかった。

このような背景から、コストと時間がかかりすぎる実機ＭＲＩ装置を必ず必要とする従来型のＭＲＩあるいはＭＲＩのパルスシークエンス開発は時代遅れの感があり、充分な性能を保持しなおかつ安価なＭＲＩシミュレータによる高効率なＭＲＩ開発が必要になってきている。

本発明の目的は、短時間でシミュレーション結果を得て、ＭＲＩパルスシーケンス開発を促進できる統合的なＭＲＩシミュレーションシステムを提供することにある。

本発明の他の目的は、核磁化の緩和時間の壁を越えて短時間でシミュレーション結果を得ることができ、実機ＭＲＩ検出部を駆動できる、統合的なＭＲＩシミュレーションシステムを提供することにある。

本発明の一態様によるＭＲＩシミュレーションシステムは、被写体を撮像するＭＲＩ検出部と、数値ファントムを撮像する仮想ＭＲＩ検出部と、ＭＲＩパルスシーケンスに基づき前記ＭＲＩ検出部を制御するＭＲＩコントローラと、前記ＭＲＩパルスシーケンスに基づき前記仮想ＭＲＩ検出部を制御するＭＲＩシミュレータと、前記ＭＲＩコントローラ又は前記ＭＲＩシミュレータから出力されるＭＲＩ信号を復号して再構成画像を生成するＭＲＩリコンストラクタと、前記ＭＲＩコントローラからのＭＲＩ信号を復号した第１の再構成画像及び／又は前記ＭＲＩシミュレータからのＭＲＩ信号を復号した第２の再構成画像を表示する表示部とを有することを特徴とする。

上述したＭＲＩシミュレーションシステムにおいて、前記ＭＲＩシミュレータは、前記ＭＲＩパルスシーケンスに基づいて単位時間毎にＭＲＩシミュレーション計算を行い、前記ＭＲＩパルスシーケンスが所定のパルスシーケンスの場合には、所定期間の間に核磁化が感じる磁場強度の変化量の積分値を演算し、演算された前記積分値に基づき、積分値を演算した前記所定時間内の前記単位時間毎のＭＲＩシミュレーション計算を割愛し、１回で前記所定期間のＭＲＩシミュレーション計算を行うようにしてもよい。

本発明によれば、被写体を撮像するＭＲＩ検出部と、数値ファントムを撮像する仮想ＭＲＩ検出部と、ＭＲＩパルスシーケンスに基づきＭＲＩ検出部を制御するＭＲＩコントローラと、ＭＲＩパルスシーケンスに基づき仮想ＭＲＩ検出部を制御するＭＲＩシミュレータと、ＭＲＩコントローラ又はＭＲＩシミュレータから出力されるＭＲＩ信号を復号して再構成画像を生成するＭＲＩリコンストラクタと、ＭＲＩコントローラからのＭＲＩ信号を復号した第１の再構成画像及び／又はＭＲＩシミュレータからのＭＲＩ信号を復号した第２の再構成画像を表示する表示部とを有するので、短時間でシミュレーション結果を得る統合的なＭＲＩシミュレーションシステムを提供することができる。

更に、本発明によれば、ＭＲＩシミュレータは、ＭＲＩパルスシーケンスに基づいて単位時間毎にＭＲＩシミュレーション計算を行い、ＭＲＩパルスシーケンスが所定のパルスシーケンスの場合には、所定期間の間に核磁化が感じる磁場強度の変化量の積分値を演算し、演算された積分値に基づき、積分値を演算した所定時間内の単位時間毎のＭＲＩシミュレーション計算を割愛し、１回で所定期間のＭＲＩシミュレーション計算を行うようにしたので、核磁化の緩和時間の壁を越えて短時間でシミュレーション結果を得ることができ、実機ＭＲＩ検出部を駆動できる、統合的なＭＲＩシミュレーションシステムを提供することができる。

本発明の一実施形態によるＭＲＩソフトウエアプラットフォームを用いたＭＲＩシミュレーションシステムを示す図である。 本発明の一実施形態によるＭＲＩシミュレーションシステムを示す図である。 本発明の一実施形態によるＭＲＩシミュレーションシステムにおける、被写体画像データを用いた静磁場の不均一性の抽出方法の説明図である。 従来のＭＲＩシミュレーションプログラムの概略を示す図である。 本発明の一実施形態によるＭＲＩシミュレーションシステムにおけるＭＲＩシミュレーションプログラムの概略を示す図（その１）である。 本発明の一実施形態によるＭＲＩシミュレーションシステムにおけるＭＲＩシミュレーションプログラムの概略を示す図（その２）である。 本発明の一実施形態によるＭＲＩシミュレーションシステムのＭＲＩシミュレータにおける計算の高速化方法の説明図である。 本発明の一実施形態によるＭＲＩシミュレーションシステムのＭＲＩシミュレータにおけるソースコードの自動生成方法のフローチャートである。 本発明の一実施形態によるＭＲＩシミュレーションシステムにおけるＭＲＩ撮像シーケンスとソースコードとの対応を示す図である。 本発明の一実施形態によるＭＲＩシミュレーションシステムの実施例１を示す図である。 本発明の一実施形態によるＭＲＩシミュレーションシステムの実施例２を示す図である。

［発明の概要］
ＭＲＩ装置（磁気共鳴映像装置）は、日進月歩での進化が宿命づけられており、相当の人数の研究開発人員と、ユーザーまで含めた、ＭＲＩ装置と撮像手法のイノベーションの枠組みを構築・維持する努力が行われている。しかしＭＲＩ装置そのものが非常に高価であり撮像に数分〜数時間の実時間が必要であることから、旧来の大規模な展開にかかわるコストは甚大で、将来に向けてはコンパクトで効率の良いＭＲＩ開発の手法が希求されている。特に、複雑なハードウェア制御を伴う「撮像パラメタの最適化」に最もコストがかかる。

本発明は、上記の問題を解決するために、ＧＰＧＰＵ（ＧＥＮＲＡＬ ＰＵＲＰＯＳＥ ＧＲＡＰＨＩＣ ＰＲＯＣＥＳＳＯＲ ＵＮＩＴ）を搭載したＰＣ（パーソナルコンピュータ）内又はネットワーク・サーバー上に仮想的なＭＲＩ検出部を実現し、ＭＲＩ実機を保有していなくても、ＭＲＩの模擬撮像を行うことができるようにするものである。

安価かつ並列化が可能なＧＰＧＰＵやクラウドによる超高速演算によって、実際にＭＲＩ装置で撮像を行うよりも短時間で、ＭＲＩシミュレータによるシミュレーション結果が得られるので、撮像パラメタの最適化の効率が飛躍的に向上する。

また、このＭＲＩシミュレータをシステム化してＭＲＩハードウェアを駆動制御できる機能を持たせることで、ＭＲＩシミュレータとＭＲＩ装置との差を最小に保つことができるようになり、将来にわたってＭＲＩ装置の開発効率を高めることができる。

これまでも、コイル単体などの個別要素の開発にシミュレータは用いることが提案されているが、本発明は、ＭＲＩ信号発生の物理原理であるＢｌｏｃｈ方程式の適用から始まりＭＲＩ撮像の全体を網羅していることに特徴がある。

また、ＭＲＩシミュレータは、実際のＭＲＩ装置に代わって、オシロスコープで観測する実際の物理現象と比較しても遜色ない疑似ＭＲＩ信号を作り出すことができるので、医学生や工学系学生の教育用途としても有用である。

［ＭＲＩシミュレーションシステム］
本発明の一実施形態によるＭＲＩシミュレーションシステムについて図１乃至図１０を用いて説明する。

（ＭＲＩソフトウエアプラットフォームを用いたＭＲＩシミュレーションシステム）
図１に本実施形態によるＭＲＩソフトウエアプラットフォームを用いたＭＲＩシミュレーションシステムを示す。

本実施形態のＭＲＩシミュレーションシステム１０は、図１に示すように、ＭＲＩソフトウエアプラットフォーム２０と、光伝送プラットフォーム３０と、実機であるＭＲＩ検出部４０により構成されている。ＭＲＩソフトウエアプラットフォーム２０と実機であるＭＲＩ検出部４０とは光伝送プラットフォーム３０を介して接続される。

ＭＲＩソフトウエアプラットフォーム２０は、ＭＲＩ装置やＭＲＩ手法の開発者のための装置的な枠組みである。ＭＲＩ 検出部４０を駆動して実際の被写体をＭＲＩ撮像することと、保存された数値化ファントムのシミュレーションを行うことができる。実際のＭＲＩ撮像によって、数値化ファントムを充実させることもできる。

光伝送プラットフォーム３０は、ＭＲＩ信号の伝送のために、光ファイバ接続を用いてＭＲＩソフトウエアプラットフォーム２０と実機であるＭＲＩ検出部４０とを接続する。信号検出器と分光器本体を光で切り分けることができるので、アナログノイズの回り込みや信号の劣化を防ぐ利点がある。

ＭＲＩソフトウエアプラットフォーム２０は、ＭＲＩシーケンスオプティマイザー２１と、ＭＲＩシーケンスジェネレータ２２と、ＭＲＩコントローラ２３と、ＭＲＩシミュレータ２４と、仮想ＭＲＩ検出部２５と、数値ファントム２６により構成されている。

ＭＲＩシーケンスオプティマイザー２１は、入力されたＭＲＩシーケンス（＊．ｓｅｑファイル）の被最適化パラメタを変化させ、撮像およびシミュレーション結果を目標値に近づける機能を有する。評価関数を有している。

ＭＲＩシーケンスジェネレータ２２は、あらかじめライブラリとして用意したＭＲＩシーケンスのひな形（＊．ｓｅｑ）を用いて、使用者から入力された撮像＆シミュレーションパラメタに適ったＭＲＩシーケンス・ファイル（＊．ｓｅｑ）を生成する機能を有する。

ＭＲＩコントローラ２３は、ＭＲＩで使用するハードウェアを実時間で制御する機能を有する。例えば、本願と同一出願による特願２０１３−５３５６９１のように、時間分解能１μ秒で各機器の制御を行う。ソフトウェア内で、ハードウェアの動きを完全に記述できるので、ＭＲＩシミュレータで使用される疑似的なハードウェアの挙動を実現するのに障壁がない。

ＭＲＩシミュレータ２４は、Ｂｌｏｃｈ方程式に従って核磁化（核スピンの集合体）を三次元的な画像マトリクスに見立て、ＭＲＩで行っている撮像すなわちＭＲＩ現象をコンピュータ内でシミュレートする機能を有する。単に、一画素内で起こるＮＭＲ現象をシミュレートする目的にも使用することができる。

仮想ＭＲＩ検出部２５は、実機であるＭＲＩ検出部４０をコンピュータ上で仮想的に実現したものである。仮想ＭＲＩ検出部２５は、ＭＲＩの検出部ハードウェアである静磁場磁石、勾配（傾斜）磁場コイル、シムコイル、送信ＲＦコイル、受信ＲＦコイルを仮想的に表現したものである。

静磁場磁石はＭＲＩ撮像を行う均一な静磁場空間を０．０５ｍＴｅｓｌａから２５Ｔｅｓｌａ程度で場合に応じて発生し、ハードウェアの不完全性による静磁場の不均一性を有している。勾配（傾斜）磁場コイルは、ＭＲＩ撮像を行うための勾配磁場空間を有しているが、中心から離れるにつれて大きくなる非線形性と、時間的に変化する渦電流による勾配磁場も有している。シムコイルは、静磁場の不均一性を打ち消して除去するために使用される電流コイルである。送信ＲＦコイルは、被写体の核磁化を励起するために用いられる。均一な領域、不均一な領域の分布がある。受信ＲＦコイルは、被写体の核磁化の信号を受信するために用いられる。均一な領域、不均一な領域の分布がある。同時に複数個用いられる場合があり、この場合、相互のカップリングが問題となる。送信ＲＦコイル、受信ＲＦコイルは、兼用されることがある。

数値ファトム２６は、仮想ＭＲＩ検出部２５で検出できるように、実際の被写体に相当するものをコンピュータ上で仮想的に実現したファントムである。単純に被写体の２次元あるいは３次元構造を画像データとして保持するだけでなく、ＮＭＲ／ＭＲＩの物理パラメタであるＴ１緩和時間、Ｔ２緩和時間、共鳴周波数、ＮＭＲ原子核の量、拡散係数等を格納した多次元データとなっている。ＮＭＲ／ＭＲＩパラメタだけでなく、リアルタイムで変化する生命現象として心拍の動きや流体現象として血流を模擬させてもよい。

数値ファントム２６には、コンピュータ上で仮想的に作成した純粋な数値ファントムと、被写体のＭＲＩデータを数値化して作成した数値化ファントムとがある。

光伝送プラットフォーム３０は、検出用アンプ３１と光伝送器３２と光伝送器３３と制御器・送受信機３４により構成されている。

検出器アンプ３１は、実機のＭＲＩ検出部４０からのアナログのＭＲＩ信号を取り込みデジタル化する。制御器・送受信機３４は、ＭＲＩコントローラ２３により制御される。

光伝送器３２が送信側の場合は、デジタル化されたＭＲＩ信号を、光伝送可能なデジタル信号に変換し、ＭＲＩプラットフォーム側の光伝送器３３に送信する。光伝送器３３が受信側の場合は、光伝送されてきたＭＲＩ信号を、制御器・送受信機３４が扱うデジタル信号の形式に変換する。制御器・送受信機３４は、ＭＲＩコントローラ２３により制御される。光伝送器３３が送信側の場合は、制御器・送受信機３４から送られてきたＭＲＩ検出部４０を駆動するための制御信号を、光伝送可能なデジタル信号に変換し、ＭＲＩ検出器側の光伝送器３２に送信する。光伝送器３２が受信側の場合は、光伝送されてきたＭＲＩ制御信号を、アンプ４２に送信する。

実機であるＭＲＩ検出部４０は、静磁場磁石（図示せず）、勾配磁場コイル（図示せず）、ＲＦコイル（図示せず）、シムコイル（図示せず）、それらを駆動する電源（アンプ４２）、付帯ユニット（図示せず）から構成される。

実機であるＭＲＩ検出部４０は、ＭＲＩ磁石４１とアンプ４２とを有している。内部にある測定対象物である被写体４３のＮＭＲ／ＭＲＩ信号を検出することができる。

ＭＲＩ検出部４０からのＮＭＲ／ＭＲＩ信号は、光伝送プラットフォーム３０の検出用アンプ３１に入力される。光伝送プラットフォーム３０の光伝送器３２からのＭＲＩ制御信号は、実機であるＭＲＩ検出部４０のアンプ４２に入力されている。

（ＭＲＩソフトウエアプラットフォーム）
図２に本実施形態によるＭＲＩシミュレーションシステムにおけるＭＲＩソフトウエアプラットフォーム２０を示す。

ＭＲＩソフトウエアプラットフォーム２０に、ＭＲＩ検出部４０および被写体４３を加えたものが、ＭＲＩシミュレーションシステムとなる。

ＭＲＩシーケンスライブラリ５０は、ＭＲＩの様々な撮像法をファイル（＊．ｓｅｑ）としてまとめたものである。教科書（Ｔｅｘｔｂｏｏｋ）や、論文（Ｐａｐｅｒ）、特許出願（Ｐａｔｅｎｔ）等の資料をもとに準備する。

ＭＲＩソフトウエアプラットフォーム２０において、ＭＲＩシーケンスジェネレータ（ＭＲＩ信号符号器）２２は、ＭＲＩシーケンスライブラリ５０からひとつのファイル（＊．ｓｅｑ）を取り出し、ＭＲＩシーケンスを生成する。ＭＲＩシーケンスは、ＭＲＩコントローラ２３又はＭＲＩシミュレータ２４を駆動させ、実際にあるいは疑似的にＭＲＩ信号を発生させる。ＭＲＩ信号とは、ＮＭＲ信号に画像情報を符号化したものである。

ＭＲＩシーケンスをＭＲＩコントローラ２３に入力すれば、実機であるＭＲＩ検出部４０により、被写体４３のＭＲＩ撮像が実施できる。

ＭＲＩシーケンスをＭＲＩシミュレータ２４に入力し、仮想ＭＲＩ検出部２５に数値ファントム２６を入力すれば装置パラメタを反映した状態でＭＲＩシミュレーションが実行できる状態となり、数値ファントム２６の仮想的ＭＲＩ撮像が実施できる。ＭＲＩシミュレータ２４は、計算の高速化のためにＧＰＧＰＵ（ｇｅｎｅｒａｌ ｐｕｒｐｏｓｅ ｇｒａｐｈｉｃ ｐｒｏｃｅｓｓｏｒ ｕｎｉｔ）を登載している。

なお、必ずしもローカルにＧＰＧＰＵハードウェアを保持している必要はなく、高速ネットワーク経由で、疑似的に実現してもよい。

ＭＲＩリコンストラクタ（ＭＲＩ信号復号器）２７は、ＭＲＩコントローラ２３又はＭＲＩシミュレータ２４の駆動により生成された各ＭＲＩ信号を、ＭＲＩシーケンスの情報から復号して画像化する、いわゆる画像再構成をして、被写体４３からは再構成画像２８Ａ、数値ファントム２６からは再構成画像２８Ｂを生成し、各画像を保存してモニタ（図示せず）に表示する。

再構成画像２８Ａは、ＭＲＩコントローラ２３からのＭＲＩ信号により再構成された実機のＭＲＩ検出部４０による撮像画像である。再構成画像２８Ｂは、ＭＲＩシミュレータ２４からのＭＲＩ信号により再構成された仮想ＭＲＩ検出部２５による撮像画像である。

再構成画像２８Ａおよび再構成画像２８Ｂは、それぞれ被写体として選択した人間やネズミのＭＲＩ画像データである。空間的には２次元ないし３次元マトリクスとして格子状に構成されており、プロトン密度強調画像、Ｔ１緩和時間強調画像、Ｔ２緩和時間強調画像等に相当する。被写体４３または数値ファントム２６に時間変化するものを選んだ場合は、画像データに時間軸が追加されることになる。

また、再構成画像２８Ａ、２８Ｂは、画像再構成処理（フーリエ変換）後のものを意図するが、画像再構成処理をせずに、いわゆる生データをそのまま表示／保存したものでもよい。

ＭＲＩコントローラ２３からのＭＲＩ信号をオシロスコープ２９Ａにより直接観測することができる。ＭＲＩシミュレータ２４からのＭＲＩ信号を仮想現実的にオシロスコープ２９Ｂにより直接観測することができる。

再構成画像２８Ａ、２８Ｂをモニタ（図示せず）に比較して表示し、ＭＲＩコントローラ２３、ＭＲＩシミュレータ２４からの各ＭＲＩ信号をオシロスコープ２９Ａ、２９Ｂにより同時に観測することにより、研究者等は、実機のＭＲＩ検出部４０の不具合や、仮想ＭＲＩ検出部２５の不具合、ＭＲＩコントローラ２３の不具合、ＭＲＩシミュレータ２４の不具合、ＭＲＩシーケンスライブラリ５０の完成度、ＭＲＩシーケンスジェネレータ２２の不具合、ＭＲＩリコンストラクタ２７の不具合、被写体４３の状態、数値ファントム２６の完成度等を把握することができる。

（実際の被写体と数値ファントム）
数値ファトム２６には、実際の被写体に相当するものをコンピュータ上で仮想的に実現した純粋な数値ファントムと、被写体のＭＲＩデータを数値化して作成した数値化ファントムとがある。

純粋な数値ファントムはコンピュータで作成する疑似的な被写体であり、人体やマウスなどの形状・構造を画素の集合体として立体的に模すだけでなく、ＭＲＩパラメタである、密度分布、水／脂の違いを示す共鳴周波数の差、Ｔ１緩和時間、Ｔ２緩和時間、拡散係数等のＭＲＩで計測できる多次元の物理値を画素に付与することで作成する。

数値ファントムには、このように単純に静的な疑似的被写体だけでなく、例えば毛細血管内を移動する血流を模して生命現象をもシミュレートした、動的な数値ファントムもある。脳梗塞やＢＯＬＤ効果等で、ＭＲＩの画像コントラストが得られる決定的な要因が解明されていないものも多く、数値ファントムを用いたＭＲＩシミュレーションを実施することが画像コントラストのメカニズム解明に大きく貢献している。

いくつかのパラメタを限定的に振ってあらかじめ得ておいたＮＭＲ／ＭＲＩシミュレーションの結果と、実際のＭＲＩ撮像の結果をマッチング比較することによって、尤もらしいパラメタを推測（計測）することができる。

数値化ファントムは、人間やマウス等の実測画像データを用いて作成する数値ファントムである。純粋な数値ファントムはシミュレーションデータとしては理想的であるが、実際に取得できる被写体のＭＲＩデータとの乖離がある。そこで、ＭＲＩ検出部で実際に取得したＭＲＩ画像を解析して、多次元化することで作成するのが数値化ファントムである。

理想的には、数値化ファントムには、プロトン密度やＴ１緩和時間、Ｔ２緩和時間の他に、計測し得るすべてのＭＲＩパラメタ（水／油、拡散係数、磁化率等）を保持していることが好ましい。数値化ファントムに、それらＭＲＩパラメタを保持することは、多くのＭＲＩシミュレーションにおいて有用である。

本願明細書では、純粋な数値ファントムと数値化ファントムを総称して数値ファントムと称している。

（数値化ファントムの利用）
図２に示すＭＲＩソフトウエアプラットフォーム２０では、被写体を撮像したＭＲＩデータを数値ファントム２６に取り入れやすくするための様々な工夫をしている。

例えば、ＭＲＩリコンストラクタ（ＭＲＩ信号復号器）２７から出力された被写体のＭＲＩデータ５１を用いて数値ファントム（数値化ファントム）２６を作成する。

また、被写体４３の実験では取得しないＭＲＩ情報５２、例えば、自己拡散係数、磁化率等を抽出し、そのＭＲＩ情報５２を用いて数値ファントム（数値化ファントム）２６を作成する。

また、後述する方法によりＭＲＩ検出部４０の静磁場の不均一性を抽出し、その静磁場の不均一性データ５３を用いて仮想検出部２５のパラメタを作成する。

このようにして、ＭＲＩ検出部４０や被写体４３の特性を反映した数値化ファントムを作成することにより、研究者等は、実際のＭＲＩ装置を長時間占有することなく、ＭＲＩ総ソフトウエアプラットフォームを用いて、ＭＲＩの装置開発あるいはＭＲＩシーケンスの開発を実行することができる。

（静磁場の不均一性の抽出方法）
図３に本実施形態によるＭＲＩシミュレーションシステムにおける、被写体画像データを用いた静磁場の不均一性の抽出方法を示す。

ＭＲＩ検出部４０における静磁場の不均一性は、おもに主磁場の超電導磁石および永久磁石の製造上避けられない誤差や磁性物によって生じるものである。理想的には０ｐｐｍであるが、実際には±１０００ｐｐｍ程度の不均一性となる場合がある。

通常、いわゆるシミングによって静磁場の均一化が行われているが、それでも残留した成分を静磁場の不均一性と呼ぶ。この静磁場の不均一性を２回のスピンエコー撮像法の位相画像の差分によって計測することができる。

スピンエコー撮像法では、一般に、ＲＦ磁場照射によるスピンエコー中心（焦点）と勾配（傾斜）磁場コイルの駆動によって発生する勾配エコー中心（焦点）を合致させた撮像が行われている。スピンエコー中心と勾配エコー中心を合致させ、スピンエコー中心と勾配エコー中心とのずれ時間であるΔｔ＝０として撮像する。これを静磁場の不均一性の計測における基準画像とする。図３の中央のＭＲＩ画像が基準画像である。

静磁場の不均一性を計測する場合は、核磁化の位相ずれを検出するために、位相差分を取得する画像のスピンエコー中心、もしくは勾配エコー中心を、基準画像の場合と比較して例えば１ｍｓ遅らせて、画像取得を行う。例えば、スピンエコー中心と勾配エコー中心との合焦時間をΔｔ＝１ｍｓずらして撮像する。実際に、撮像シーケンスを調整する場合には勾配エコーが見えている状態ではスピンエコーは目視できないため、勾配磁場の印加をＯＮ／ＯＦＦした状態で、スピンエコー中心と勾配エコー中心を識別する必要がある。図３の左側のＭＲＩ画像が取得した画像である。静磁場の不均一性を反映している。

これらのＭＲＩ画像を、核磁化の方位を示す位相画像に変換し、必要に応じて位相アンラッピング（ｕｎｗｒａｐｐｉｎｇ）処理を施し、位相差画像を生成することで静磁場の不均一性を抜き出すことができる。図３の右側のＭＲＩ画像が静磁場の不均一性のみが抽出された画像である。

位相差画像から不均一性への計算式は下記のようになる。例えば、計算された位相差画像内のある点において核磁化の位相差が３６０度であった場合、Δｔ＝１ｍｓを選んでいた場合は、核磁化の１ｍｓに回転した量に相当するので、１秒に換算すると、１０００回転つまり１ｋＨｚの静磁場の不均一性が存在していることになる。なお、撮像時に設定する遅延のΔｔ＝１ｍｓは任意であり、計測する静磁場の不均一性と静磁場強度を勘案して、通常０．１ｍｓから２ｍｓ程度の範囲で選択される。

（従来のＭＲＩシミュレーションプログラム）
図４に従来のＭＲＩシミュレーションプログラムの概略を示す。

ＭＲＩシミュレーションプログラム６０には、三次元的な被写体の特性情報のパラメタである、密度分布、Ｔ１マップ（Ｔ１緩和時間分布）、Ｔ２マップ（Ｔ２緩和時間分布）、静磁場分布（中心共鳴周波数からのシフト分の分布）のパラメタ６１と、ＭＲＩ検出部 ４０の装置のパラメタである、勾配磁場強度／分布、静磁場強度／分布、高周波磁場強度／分布、渦電流の特性等のシステムパラメタ６２と、ＭＲＩシーケンスである、パルスシーケンス６３とが入力される。例えば、パラメタ６１の静磁場強度／分布は、静磁場の不均一性５３を反映したものである。あらかじめ、被写体の形状や磁化率によって引き起こされる静磁場の不均一性を加味して反映させてもよい。

従来のＭＲＩシミュレーションプログラム６０では、位相エンコードループ６０Ａ内にシーケンスイベントループ６０Ｂが包含されている。

シーケンスイベントループ６０Ｂでは、シーケンスのイベントを解釈し、分岐先を決定し、イベントの時刻と引数を分岐先の関数へ受け渡すステップ６０Ｂ１と、核磁化の仰角を三次元方向に変化させ、核磁化強度を減衰させるステップ６０Ｂ２と、核磁化の位相を二次元面内で変化させ、核磁化強度を減衰させるステップ６０Ｂ３と、核磁化成分の総和によりＭＲＩ信号を計算するステップ６０Ｂ４とを繰り返し実行する。

これらを実行することにより疑似的なＭＲＩ信号が出力される。

このように従来のＭＲＩシミュレーションプログラムでは、時間的に連続したシーケンスのイベントループで単位時間刻みごとにイベントを逐次的に解釈して、対応するＣＰＵ／ＧＰＧＰＵの並列計算関数を実行している。ＣＰＵ／ＧＰＧＰＵの並列計算関数は、内部で核磁化ベクトルの値、Ｔ１値、Ｔ２値等をグローバルメモリから読み書きする必要があり、データ転送の遅延待ちにより、計算処理が律速される。このため、シミュレーションに要する時間が膨張する。具体的には次のような理由である。

従来のＭＲＩシミュレーションプログラムでは、ＭＲＩ装置が実際に信号取得を行う場合の、静磁場中に静置された被写体数値化データの核磁化の挙動を、コンピュータ内で再現する。ＭＲＩは核磁化分布を可視化する画像計測なので、被写体を、例えば、１０μｍ〜１ｍｍ程度大きさの２５６^３個の隣接する直方体マトリクスとして分解して、ＮＭＲ現象の標準であるＢｌｏｃｈ方程式に従って、１μ秒から数１０μ秒の刻みで各画素のＮＭＲ現象をシミュレートする。被写体の核磁化分布（２５６^３個）と、実際に取得できるＭＲＩデータ（２５６^３画素）には、フーリエ変換の関係がある。このためＭＲＩデータを２５６^３画素で収集する場合は、２５６^３個のＮＭＲ現象シミュレーションを２５６^３回以上繰り返す必要があり、膨大な処理時間を必要とする。

なお、各画素のＮＭＲパラメタである、プロトン密度分布、Ｔ１分布、Ｔ２分布、共鳴周波数分布（不均一磁場分布）が、画像コントラストだけでなく画像の性質を決定する要素となる。

（本実施形態のＭＲＩシミュレーションプログラム）
図５及び図６に本実施形態のＭＲＩシミュレーションプログラムの概略を示す。

ＭＲＩシミュレーションプログラム７０には、従来のＭＲＩシミュレーションプログラムと同様に、三次元的な被写体の特性情報のパラメタである、密度分布、Ｔ１マップ（Ｔ１緩和時間分布）、Ｔ２マップ（Ｔ２緩和時間分布）、静磁場分布（中心共鳴周波数からのシフト分の分布）のパラメタ６１と、ＭＲＩ検出部のパラメタである、勾配磁場強度／分布、高周波磁場強度／分布、渦電流の特性等のシステムパラメタ６２と、ＭＲＩシーケンスである、パルスシーケンス６３とが入力される。

ＭＲＩシミュレーションプログラム７０では、まず、パルスシーケンス６３とシステムパラメタ６２に依存したＣＵＤＡ（クーダ）のカーネル関数のソースコードを自動生成する（ステップ７０Ａ）。外部から入力したＭＲＩパルスシーケンスの動作であるイベントのコードを解釈し、外部から入力したＭＲＩシステムのパラメタを用いて、ＧＰＧＰＵの性能を活かした実行コードとして実装されているカーネル関数のソースコードを、パルスシーケンスのプログラムが最適化されるようなアルゴリズムで自動的に文字列として生成する。

ＣＵＤＡ（Ｃｏｍｐｕｔｅ Ｕｎｉｆｉｅｄ Ｄｅｖｉｃｅ Ａｒｃｈｉｔｅｃｔｕｒｅ：クーダ）とは、ＮＶＩＤＩＡ Ｃｏｒｐｏｒａｔｉｏｎ（エヌビディアコーポレーション）により提供されるＧＰＧＰＵ向けのＣ言語の統合開発環境であり、コンパイラ（ｎｖｃｃ）やライブラリなどから構成されている。アプリケーションを実行する基盤となるプラットフォーム／アーキテクチャそのものをＣＵＤＡと呼ぶこともある。

次に、ＮＶＩＤＩＡ Ｃｏｒｐｏｒａｔｉｏｎ（エヌビディアコーポレーション）により提供されるランタイムコンパイルライブラリによりＣＵＤＡのカーネル関数のソースコードをコンパイルしてコンピュータにより実行可能なオブジェクトコード（機械語）に変換する（ステップ７０Ｂ）。例えば、ＧＰＧＰＵの製造メーカが提供する、ＧＰＧＰＵのソースコードのコンパイラを用いて、ＭＲＩシミュレーションプログラム上から、ソースコードを実行可能コードに変換する。

次に、ＣＵＤＡカーネル関数によるブロッホ（Ｂｌｏｃｈ）方程式の計算とＭＲＩ信号の収集を行う（ステップ７０Ｃ）。ブロッホ（Ｂｌｏｃｈ）方程式とは、磁気共鳴の現象論的記述をする方程式である。

ステップ７０Ｃは、実際のプログラム実行プロセスであり、ＧＰＧＰＵのホストとなるＣＰＵとＧＰＧＰＵが同期をとりながら、核磁化シミュレーションを高速に実行する。

ステップ７０Ｃは、図６に示すように、位相エンコードループ７０Ｃ１内に、実行時にコンパイルされるＣＵＤＡ関数７０Ｃ２が内包され、関数７０Ｃ２内に、ＧＰＧＰＵで並列化された画素計算ループ７０Ｃ３が内包されている。

位相エンコードループ７０Ｃ１は、位相エンコードループであり、核磁化の位置情報を付加するために、核磁化シミュレーションで行われる最も外側のループである。三次元イメージングでは、画素数に対応したループであり、通常、二重のループとなっていて、本実施形態では、ＣＰＵで実行される。

ＣＵＤＡ関数７０Ｃ２は、核磁化の運動を逐一計算する、ＧＰＧＰＵの演算コード（カーネル関数）が実行順に並べられている。

画素計算ループ７０Ｃ３には、Ｒｏｔａｔｅ３ＤａｎｄＲｅｌａｘ（…）；Ｒｏｔａｔｅ２ＤａｎｄＲｅｌａｘ（…）；ＡｃｑｕｉｒｅＮＭＲＳｉｇｎａｌ（…）等が含まれている。

画素計算ループ７０Ｃ３は、カーネル関数を核磁化計算の機能毎に並べたものである。Ｒｏｔａｔｅ３ＤａｎｄＲｅｌａｘ（）という関数は、不均一磁場中における高周波パルスによる励起、すなわち核磁化が緩和しながら三次元的に回転する操作を表す。Ｒｏｔａｔｅ２ＤａｎｄＲｅｌａｘ（）は、不均一磁場中で、高周波パルスが印加されていないときに、緩和しながらＺ軸の周りに歳差運動する操作を表す。ＡｃｑｕｉｒｅＭＲＩＳｉｇｎａｌ（）は、歳差運動している核磁化の成分をすべて足し合わせることにより、ＭＲＩ信号を計算する操作を表す。

これらを実行することにより擬似的なＭＲＩ信号が出力される。

このように本実施形態のＭＲＩシミュレーションプログラムでは、シーケンスとシステムパラメタが入力された時点で、計算手順は決定する。従来のＭＲＩシミュレーションプログラムでは逐次的に分岐判断しながら行っていた計算を、分岐のないソースコードとして生成する。それをコンパイルすることにより、従来のＭＲＩシミュレーションプログラムのシーケンスイベントループから分岐を除去することができる。また、それぞれの計算は、レジスタ内で行うことができ、ループの先頭と末尾でメモリ読み書きをするだけでよいため、メモリ読み書きによる性能低下も防げる。このため、シミュレーションに要する時間が大幅に短縮される。

（ＭＲＩシミュレーション計算の高速化方法）
図７に本実施形態のＭＲＩシミュレーションシステムのＭＲＩシミュレータにおける計算の高速化方法を示す。この計算の高速化方法は、特に、勾配磁場の印加と、静磁場の不均一性によって、核磁化ベクトルが変化する現象に関係する。

本実施形態のシミュレーション計算の高速化方法は、ＭＲＩのパルスシークエンスを時間軸上でＭＲＩ制御イベントの特徴毎に分割し、「勾配磁場と静磁場の不均一性に関する計算」を行う際には、Δ磁場×Δ時間の積分値を適用することで連続した単位時間当たりのシミュレーションを不要として計算を高速化している。

通常、時間変化するＭＲＩ現象をシミュレーションする場合、最小単位時間を設定してその刻みで対象の計算を行う必要がある。ＭＲＩシミュレーションにおいても同様である。

ＭＲＩシミュレーションで扱う磁場を、静磁場系をＢ０、高周波磁場系をＢ１とするならば、静磁場の不均一性はΔＢ０、勾配磁場の変化はΔＧＸ／ΔＧＹ／ΔＧＺである。静磁場系の変化量としてΔＡ＝ΔＢ０＋ΔＧＸ／ΔＧＹ／ΔＧＺ＋渦電流効果と表記することもできる。

ＭＲＩシミュレーションにおいて、場合によっては、静磁場強度Ｂ０が２３．４Ｔｅｓｌａもありうるが、１ＧＨｚに近い周波数で歳差運動する核磁化をＢ０＋ΔＡの強力な磁場で数値計算するのは現在の計算機では不可能である。

そこで本発明においてはＢ０を一定としてΔＡのみをとりあつかう回転座標系で起こるＮＭＲおよびＭＲＩ現象のシミュレーションを用いている。

本実施形態では、シミュレーション対象であるＭＲＩパルスシーケンスの特徴に着目して単位時間ごとの計算を一気に割愛する計算方法を用いている。

上述のΔＡの勾配磁場、渦電流、および静磁場の不均一性は空間的に変化する特徴があるが、微視的な１画素においては、勾配磁場と渦電流は時間変化するが、均一磁場、静磁場の不均一性は時間変化しない均一磁場とみなすことができる。どちらの場合であっても、磁場の変化量の積分値が時間軸で得られていれば、それを用いることで、繰り返し計算を行う必要はなく、１回の計算で、高速化できる時間幅を一気に短縮することができる。磁場の変化量の積分値の計算は、ＭＲＩのシミュレーションの１ステップと比較して、非常に高速に計算結果を得ることができるので、この手法による高速化のメリットは大きい。

図７（ａ）のように核磁化が感じる磁場の強度Ａの時間変化がない場合、従来のＭＲＩシミュレーションでは、最小単位時間刻み毎にＮ回の数値計算を行っていたが、本実施形態のＭＲＩシミュレーションでは、核磁化が感じる磁場の強度Ａが一定として、その面積ΔＡ×ｔに基づいて１回の数値計算を行うだけでよい。これにより数値計算の高速化を図ることができる。

図７（ｂ）のように核磁化が感じる磁場の強度Ａの時間変化がある場合、従来のＭＲＩシミュレーションでは、最小単位時間刻み毎にＮ回の数値計算を行っていたが、本実施形態のＭＲＩシミュレーションでは、核磁化が感じる磁場の強度Ａが変化しても、その平均値に基づく積分値（面積）ΔＡ×ｔに基づいて１回の数値計算を行うだけでよい。これにより数値計算の高速化を図ることができる。

なお、図７では、時間変化がない場合として時間軸の計算を大きく割愛する例を記載しているが、ΔＡ×ｔ＝Ａ×ΔｔとなるＡを求めて、ステップ状に変化するように配置して残りは強度Ａを０として計算を割愛してもよい。

（ソースコードの自動生成）
図８のフローチャートを用いて、本実施形態のＭＲＩシミュレーションプログラム７０におけるステップ７０Ａの詳細を説明する。

ここでは特に図７で説明した高速化手法を、ＭＲシミュレータに入力されたＭＲＩシーケンスの変形に利用している。

ステップ７０Ａでは、パルスシーケンス６３とシステムパラメタ６２に依存したＣＵＤＡのカーネル関数のソースコードを自動生成する。

まず、パルスシーケンス６３から、今回ソースコードを生成する一連のパルスシーケンスを読み込む（ステップＳ１）。

次に、核磁化ベクトルの位相の積分値、直前のＧＸ，ＧＹ，ＧＺのイベント、ソースコード生成時刻を初期化する（ステップＳ２）。

次に、一連のパルスシーケンスの全てのイベントについてのループのステップＳ３を実行する。ステップＳ３は、次に説明するステップＳ４〜Ｓ９によって構成されている。

まず、ステップＳ４ではイベントの種類を判断し、その判断結果に基づいて、ステップＳ５〜Ｓ９を実行する。

画像マトリクス内の各点で表現される核磁化ベクトルは、大きさ、仰角、および位相情報を保持しており、以降、計測開始から現在の履歴をすべて表現する意を込めて静磁場に直交する回転座標内において、核磁化ベクトルの方位を、核磁化位相の積分値と表現する。具体的に採り得る値は、±１８０度であるが、１８０度を超えて、例えば、１８０度＋３６０度を５４０度と表現し保持してもよい。

ステップＳ５では、核磁化位相の積分値（計測開始から現在時刻まで）に，直前のＧＸイベントのイベントからここまでのＧＸイベントの変化量の積分値を足し合わせる。

図７（ｂ）右側で示した例を推し進め、この変化量を計算して求めたＧＸイベントを新たなＧＸのイベントとして対応する時刻に１つ配置する。直前のＧＸイベントのイベントからここまでのＧＸイベントのＧＸの変化量は１つのイベントの集約できる。集約されて排除した計算部分は実行しないので計算を高速化できる。

ステップＳ６では、ステップＳ５をＧＹについて行う。

ステップＳ７では、ステップＳ５をＧＺについて行う。

ステップＳ８では、核磁化ベクトルの位相の積分値に直前のソースコード生成時刻からこのイベントの時刻までの不均一磁場（静磁場の不均一性）の積分値を足し合わせて、２次元回転のソースコードを生成し，核磁化ベクトルの位相の積分値を０にする直前のソースコード生成時刻と、このイベントの時刻までの緩和のソースコードを生成直前のＧＸ、ＧＹ、ＧＺのイベントの時刻、ソースコード生成時刻をこのイベントの時刻に設定３次元の回転のソースコードを生成する。

ステップＳ９では、核磁化ベクトルの位相の積分値に直前のソースコード生成時刻からこのイベントの時刻までの不均一磁場の積分値を足し合わせて、位相の２次元回転のソースコードを生成し、位相の積分値を０にする直前のソースコード生成時刻と、このイベントの時刻までの緩和のソースコードを生成直前のＧＸ、ＧＹ、ＧＺのイベントの時刻、ソースコード生成時刻をこのイベントの時刻に設定したＭＲＩ信号サンプリングのソースコードを生成する。

一連のパルスシーケンスの全てのイベントについてのループが終了する（ステップＳ１０）と、引き続いて、核磁化ベクトルの位相の積分値に直前のソースコード生成時刻からＴＲまでの不均一磁場（静磁場の不均一性）の積分値を足し合わせて、位相の２次元回転のソースコードを生成直前のソースコード生成時刻からＴＲまでの核磁化の緩和のソースコードを生成する（ステップＳ１１）。

このようにして入力されたパルスシーケンスに基づいたソースコードを自動生成するステップＳ７０Ａを終了する。

（シーケンスとソースコードの具体例）
図９に本実施形態のＭＲＩシミュレーションプログラムを用いてシーケンスからソースコードを生成した具体例を示す。

図９左部に示すＭＲＩパルスシーケンスについて説明する。

これは繰り返し時間ＴＲ内で絶対時間表記されたものであり、左から実行時刻、イベント、引数に対応する。例えば、００．０１０．０００．０は１０．０ｍｓを起点にＲＦパルスの０００２番（組込みＳＩＮＣ波形の１ｍｓ幅）を照射せよ、というイベントに相当する。

このＭＲＩパルスシーケンスは、スピンエコー法のものである。このＭＲＩパルスシーケンスのＢｌｏｃｈシミュレーションを計算する場合、枠で囲った７個の特徴を自動あるいは手動で認識し、計算ステップ(1)〜(7)に分割して計算するソースコードを生成する。これをあらかじめ決めた回数であるＮ回分を繰り返し連続的に実行する。繰り返し時間ＴＲの設定が充分に長くない場合は、核磁化ベクトルの位相や大きさの履歴を引きずりながらのＢｌｏｃｈシミュレーションとなる。

計算ステップ(1)、(3)、(5)、(6)、(7)（所定のパルスシーケンス）は、核磁化ベクトルに、その区間のＧＸ、ＧＹ、ＧＺによる勾配磁場と、不均一磁場（静磁場の不均一性）の各積分値に応じたＺ方向を軸とした２次元の回転行列とその区間のＴ１、Ｔ２緩和行列を乗じ、Ｔ１緩和ベクトルを加算するソースコードを生成する計算ステップである。

計算ステップ(2)、(4)は、核磁化ベクトルに、その区間のＧＸ、ＧＹ、ＧＺによる勾配磁場と、不均一磁場（静磁場の不均一性）、これらに直交するＲＦ高周波磁場（Ｂ１）の積分値に応じた、核磁化ベクトルの方位に関する３次元の回転行列とその区間のＴ１、Ｔ２緩和行列を乗じ、Ｔ１緩和ベクトルを加算するソースコードを生成する計算ステップである。計算精度を確保するために３次元の回転と緩和を何分割化して計算するソースコードを生成する。

計算ステップ(6)′（所定のシミュレーション計算）は、ＭＲＩにおけるＡＤ変換器によるデータサンプリングとその前後時間に相当する。核磁化ベクトルの回転と緩和の計算をしながら、ＮＭＲ信号として、核磁化を足し合わせる計算をするソースコードを生成する計算ステップである。

この例のように、ＡＤサンプリング中の場合には、計算ステップ(6)（所定のパルスシーケンス）から分岐した特別に高速なソースコード(6)′（所定のシミュレーション計算）を生成する。計算ステップ(6)′におけるＮＭＲＳｉｇｎａｌＳａｍｐｌｉｎｇ（）関数による核磁化ベクトルの超高速な計算結果は、繰り返し時間ＴＲ内の決められた区間をサンプリングして、ＭＲＩ画像の生データの一部として別途保存した後、次の計算ステップ(7)のシミュレーションデータとしては継承しないで捨てる。

ＮＭＲＳｉｇｎａｌＳａｍｐｌｉｎｇ（）関数が最小時間分解能で行う連続した超高速な数値計算は誤差分を多く含むので、計算ステップ(6)は並行して上述の計算ステップ(1)、(3)、(5)、(7)と同様の計算を行っておき、あるいは、後の計算ステップ(7)にて、(6)と(7)の勾配磁場変化をまとめて計算するソースコードを生成する。このように切り捨てる分岐経路を持つことで、全体的な高速化と部分的な高速化を、計算精度を保ちながら両立させることができる。

図９右部に示すソースコード記述（呼出関数）について説明する。

計算ステップ(1)の呼出関数は、最初のＲＦパルスまでの不均一磁場による核磁化の２次元回転と緩和の計算を行う呼出関数である。

計算ステップ(2)の呼出関数は、ＲＦ磁場と不均一磁場による核磁化の３次元回転と緩和の計算をＮ分割して行う呼出関数である。

計算ステップ(3)の呼出関数は、ＧＸの勾配磁場と次のＲＦパルスまでの不均一磁場による核磁化の２次元回転と緩和の計算を行う呼出関数である。

計算ステップ(4)の呼出関数は、ＲＦ磁場と不均一磁場による核磁化の３次元回転と緩和の計算をＮ分割して行う呼出関数である。

計算ステップ(5)の呼出関数は、ＧＹ，ＧＺの位相エンコード勾配磁場とサンプリングをするまでのリード勾配磁場と不均一磁場による核磁化の２次元回転と緩和の計算を行う呼出関数である。

計算ステップ(6)の呼出関数は、ＧＸの勾配磁場と次のＲＦパルスまでの不均一磁場による核磁化の２次元回転と緩和の計算を行う呼出関数である。

計算ステップ(6)′の呼出関数は、リード勾配磁場と不均一磁場による２次元の回転と緩和の計算を行い、ＮＭＲ信号の計算を行う呼出関数である。信号を計算した後の核磁化は、保存せずに捨てる。

計算ステップ(7)の呼出関数は、サンプリング中のリード勾配磁場と残りのＴＲまでの不均一磁場による２次元回転と緩和の計算を行う呼出関数である。

図９に示すＭＲＩシミュレーションの具体例に適用する場合について説明する。このＭＲＩシミュレーションでは、枠で囲った７個の計算ステップ(1)〜(7)に分割して計算するソースコードを生成する。

７個の計算ステップ(1)〜(7)のうち、計算ステップ(1)、(3)、(5)、(6)、(7)は、ＲＦ高周波磁場による影響がない場合の、核磁化ベクトルの位相のシミュレーション計算であるので、本実施形態による計算の高速化方法を適用できる。

ただし、ＡＤサンプリングに相当する計算ステップ(6)については、並行して高精細なシミュレーションを切り捨てで実施する。

計算ステップ(1)、(3)、(5)、(6)、(7)は、核磁化ベクトルに、その区間のＧＸ、ＧＹ、ＧＺによる磁場と付帯する渦電流、不均一磁場（静磁場の不均一性）、の積分値に応じたＺ方向を軸とした２次元の回転行列とその区間のＴ１、Ｔ２緩和行列を乗じ、Ｔ１緩和ベクトルを加算するソースコードを生成する計算ステップである。

従来のＭＲＩシミュレーションが、ΔＡ＝勾配磁場＋渦電流＋不均一磁場について最小時間単位で起点から終点までのすべての時刻において数値計算していた「全道程なぞり型」とするならば、本実施形態では、起点から終点を一直線につないだ「直線移動型」あるいは、起点から終点までの変化量を０として瞬間移動し、終点にすべての変化量（積分値）を上乗せして帳尻を合わせた「瞬間移動型」として考えることができる。

図９のシーケンスチャートにおいて、勾配磁場のステップ変化のみを加味する場合、図７（ａ）の高速化手法を計算ステップ(1)、(3)、(5)、(6)、(7)に適用することができる。

図９のシーケンスチャートにおいて、勾配磁場の緩やかな変化、渦電流、不均一磁場の影響を加味する場合、図７（ｂ）の高速化手法を計算ステップ(1)、(3)、(5)、(6)、(7)に適用することができる。

ただし、計算ステップ(6)に関しては、ＡＤサンプリングのための詳細な計算が必要であるため、計算ステップ(6)の直前まで繋がれてきた計算結果を用いて、繰り返し時間ＴＲ毎にＡＤサンプリングのための高速核磁化計算(6)′を行う。このため、図９のソースコードの呼出関数の記載において「ＮＭＲＳｉｇｎａｌＳａｍｐｌｉｎｇ（）；」という特別な関数を用意している。

（実施例）
本実施形態のＭＲＩシミュレーションシステムの実施例を図１０及び図１１に示す。

図１０は、球状の空孔を持つファントム（模擬物体）のシミュレーションと実験画像とを比較した実施例１である。実験画像に対し，３２７６．８／２９８．６＝約１１倍の速度で、誤差のほとんどないシミュレーション画像が得られた。

図１１は、円筒状の空孔を持つファントム（模擬物体）のシミュレーションと実験画像とを比較した実施例２である。実験画像に対し、８１９．２／１７．３＝約４７倍の速度で、誤差のほとんどないシミュレーション画像が得られた。

［変形実施形態］
本発明は上記実施形態に限らず種々の変形が可能である。

例えば、上記実施形態では、ＭＲＩシミュレータを単一としたが、複数のＭＲＩシミュレータを設けてもよい。

例えば、上記実施形態では、ＭＲＩソフトウエアプラットフォームを単一としたが、複数のＭＲＩソフトウエアプラットフォームを配置／接続するようにしてもよい。

１０ ＭＲＩシミュレーションシステム
２０ ＭＲＩソフトウエアプラットフォーム
２１ ＭＲＩシーケンスオプティマイザー
２２ ＭＲＩシーケンスジェネレータ（ＭＲＩ信号符号器）
２３ ＭＲＩコントローラ
２４ ＭＲＩシミュレータ
２５ 仮想ＭＲＩ検出部
２６ 数値ファントム
２７ ＭＲＩリコンストラクタ（ＭＲＩ信号復号器）
２８Ａ、２８Ｂ 再構成画像
２９Ａ、２９Ｂ オシロスコープ
３０ 光伝送プラットフォーム
３１ 検出用アンプ
３２ 光伝送器
３３ 光伝送器
３４ 制御器・送受信機
４０ ＭＲＩ検出部
４１ ＭＲＩ磁石
４２ アンプ
４３ 被写体
５０ ＭＲＩシーケンスライブラリ
５１ 被写体のＭＲＩデータ
５２ 被写体のＭＲＩ情報
５３ 静磁場の不均一性
６０ ＭＲＩシミュレーションプログラム
６０Ａ 位相エンコードループ
６０Ｂ シーケンスイベントループ
６１ 被写体の特性情報のパラメタ
６２ ＭＲＩ検出部のパラメタ
６３ ＭＲＩパルスシーケンス
７０ ＭＲＩシミュレーションプログラム

Claims (6)

1. 被写体を撮像するＭＲＩ検出部と、
   数値ファントムを撮像する仮想ＭＲＩ検出部と、
   ＭＲＩパルスシーケンスに基づき前記ＭＲＩ検出部を制御し、ＭＲＩ信号を出力するＭＲＩコントローラと、
   前記ＭＲＩパルスシーケンスに基づき前記仮想ＭＲＩ検出部を制御し、疑似ＭＲＩ信号を出力するＭＲＩシミュレータと、
   前記ＭＲＩコントローラから出力される前記ＭＲＩ信号を復号して第１の再構成画像を生成し、前記ＭＲＩシミュレータから出力される前記疑似ＭＲＩ信号を復号して第２の再構成画像を生成するＭＲＩリコンストラクタと、
   前記ＭＲＩコントローラからの前記ＭＲＩ信号と前記ＭＲＩシミュレータからの前記疑似ＭＲＩ信号を表示し、前記ＭＲＩコントローラからの前記ＭＲＩ信号を復号した第１の再構成画像と前記ＭＲＩシミュレータからの前記疑似ＭＲＩ信号を復号した第２の再構成画像を比較して表示する表示部であって、前記ＭＲＩコントローラからの前記ＭＲＩ信号と前記ＭＲＩシミュレータからの前記疑似ＭＲＩ信号とを、前記第１の再構成画像と前記第２の再構成画像と同時観測できるように表示する表示部と
   を有することを特徴とするＭＲＩシミュレーションシステム。
2. 請求項１記載のＭＲＩシミュレーションシステムにおいて、
   前記ＭＲＩコントローラから出力される前記被写体の２次元あるいは３次元構造のＭＲＩ信号に基づいて前記被写体の２次元あるいは３次元構造の数値ファントムを作成する
   ことを特徴とするＭＲＩシミュレーションシステム。
3. 請求項１又は２記載のＭＲＩシミュレーションシステムにおいて、
   前記ＭＲＩシミュレータは、
   前記ＭＲＩパルスシーケンスに基づいて単位時間毎にＭＲＩシミュレーション計算を行い、
   前記ＭＲＩパルスシーケンスが所定のパルスシーケンスの場合には、所定期間の間に核磁化が感じる磁場強度の変化量の積分値を演算し、演算された前記積分値に基づき、積分値を演算した前記所定時間内の前記単位時間毎のＭＲＩシミュレーション計算を割愛し、１回で前記所定期間のＭＲＩシミュレーション計算を行う
   ことを特徴とするＭＲＩシミュレーションシステム。
4. 請求項３記載のＭＲＩシミュレーションシステムにおいて、
   前記ＭＲＩシミュレータは、
   前記ＭＲＩパルスシーケンスが前記所定のパルスシーケンスであって、最小時間分解能で実行することが望ましい所定のシミュレーション計算を含む場合には、前記所定のシミュレーション計算のための分岐を設け、その分岐において最小時間分解能で前記所定のシミュレーション計算を実行し、その計算結果を前記所定のパルスシーケンスのシミュレーションデータとして使用する
   ことを特徴とするＭＲＩシミュレーションシステム。
5. 請求項１乃至４のいずれか１項に記載のＭＲＩシミュレーションシステムにおいて、
   前記ＭＲＩシミュレータは、
   前記ＭＲＩパルスシーケンスに基づくカーネル関数のソースコードを生成する手段と、
   前記カーネル関数のソースコードをコンパイルしてオブジェクトコードに変換する手段と、
   前記オブジェクトコードを実行して疑似ＭＲＩ信号を出力する手段とを有する
   ことを特徴とするＭＲＩシミュレーションシステム。
6. 請求項１乃至５のいずれか１項に記載のＭＲＩシミュレーションシステムで用いられるプログラムであって、
   コンピュータを、
   前記ＭＲＩパルスシーケンスに基づくカーネル関数のソースコードを生成する手段と、
   前記カーネル関数のソースコードをコンパイルしてオブジェクトコードに変換する手段と、
   前記オブジェクトコードを実行して疑似ＭＲＩ信号を出力する手段と
   を有するＭＲＩシミュレーション装置
   として機能させるためのプログラム。