知財求人 - 知財ポータルサイト「IP Force」
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2025004325
(43)【公開日】2025-01-15
(54)【発明の名称】磁気共鳴数値シミュレーション装置、方法及びプログラム
(51)【国際特許分類】
A61B 5/055 20060101AFI20250107BHJP
【ＦＩ】
A61B5/055 311
【審査請求】未請求
【請求項の数】17
【出願形態】ＯＬ
(21)【出願番号】P 2023103933
(22)【出願日】2023-06-26
(71)【出願人】
【識別番号】594164542
【氏名又は名称】キヤノンメディカルシステムズ株式会社
(74)【代理人】
【識別番号】110003708
【氏名又は名称】弁理士法人鈴榮特許綜合事務所
(72)【発明者】
【氏名】竹島 秀則
【テーマコード（参考）】
4C096
【Ｆターム（参考）】
4C096AB26
4C096DB06
4C096DB20
(57)【要約】
【課題】磁気共鳴数値シミュレーションの計算効率を向上すること。
【解決手段】 実施形態に係る磁気共鳴数値シミュレーション装置は、入力部と数値計算部とを有する。入力部は、アイソクロマットについて、磁化の更新前の数値と磁化に対する空間方向及び／又は角周波数方向各々での偏微分の更新前の数値とを入力する。数値計算部は、磁化の更新前の数値と偏微分の更新前の数値とのうちの一部の数値を使用して、磁気共鳴の時間発展を表す更新式を演算し、アイソクロマットについて、磁化の更新後の数値と磁化に対する空間方向及び／又は角周波数方向のうちの算出対象方向での空間偏微分の更新後の数値とを算出する。
【選択図】 図１
【解決手段】 実施形態に係る磁気共鳴数値シミュレーション装置は、入力部と数値計算部とを有する。入力部は、アイソクロマットについて、磁化の更新前の数値と磁化に対する空間方向及び／又は角周波数方向各々での偏微分の更新前の数値とを入力する。数値計算部は、磁化の更新前の数値と偏微分の更新前の数値とのうちの一部の数値を使用して、磁気共鳴の時間発展を表す更新式を演算し、アイソクロマットについて、磁化の更新後の数値と磁化に対する空間方向及び／又は角周波数方向のうちの算出対象方向での空間偏微分の更新後の数値とを算出する。
【選択図】 図１
【特許請求の範囲】
【請求項1】
アイソクロマットについて、磁化の更新前の数値と前記磁化に対する空間方向及び／又は角周波数方向各々での偏微分の更新前の数値とを入力する入力部と、
前記磁化の更新前の数値と前記偏微分の更新前の数値とのうち一部の数値を使用して、磁気共鳴の時間発展を表す更新式を演算し、前記アイソクロマットについて、前記磁化の更新後の数値と前記磁化に対する前記空間方向及び／又は前記角周波数方向のうちの算出対象方向での前記偏微分の更新後の数値とを算出する数値計算部と、
を具備する磁気共鳴数値シミュレーション装置。
前記磁化の更新前の数値と前記偏微分の更新前の数値とのうち一部の数値を使用して、磁気共鳴の時間発展を表す更新式を演算し、前記アイソクロマットについて、前記磁化の更新後の数値と前記磁化に対する前記空間方向及び／又は前記角周波数方向のうちの算出対象方向での前記偏微分の更新後の数値とを算出する数値計算部と、
を具備する磁気共鳴数値シミュレーション装置。
【請求項2】
前記数値計算部は、
前記偏微分の更新前の数値を使用せず、前記磁化の更新前の数値を使用して、前記更新式を演算し、前記磁化の更新後の数値を算出し、
前記空間方向及び／又は前記角周波数方向のうちの非算出対象方向での前記偏微分の更新前の数値を使用せず、前記算出対象方向での前記偏微分の更新前の数値を使用して、前記更新式を演算し、前記算出対象方向での前記偏微分の更新後の数値を算出する、
請求項１記載の磁気共鳴数値シミュレーション装置。
前記偏微分の更新前の数値を使用せず、前記磁化の更新前の数値を使用して、前記更新式を演算し、前記磁化の更新後の数値を算出し、
前記空間方向及び／又は前記角周波数方向のうちの非算出対象方向での前記偏微分の更新前の数値を使用せず、前記算出対象方向での前記偏微分の更新前の数値を使用して、前記更新式を演算し、前記算出対象方向での前記偏微分の更新後の数値を算出する、
請求項１記載の磁気共鳴数値シミュレーション装置。
【請求項3】
前記更新式は、互いに分離された磁化更新式と空間偏微分更新式と角周波数偏微分更新式とを有し、
前記数値計算部は、
前記磁化の更新後の数値を算出するために、前記磁化の更新前の数値を参照して前記磁化更新式を演算し、
前記算出対象方向での前記空間方向の前記偏微分の更新後の数値を算出するために、前記磁化の更新前の数値及び前記算出対象方向での前記空間方向の前記偏微分の更新前の数値を参照して前記空間偏微分更新式を演算し、
前記角周波数方向の前記偏微分の更新後の数値を算出するために、前記磁化の更新前の数値及び前記角周波数方向の前記偏微分の更新前の数値を参照して前記角周波数偏微分更新式を演算する、
請求項２記載の磁気共鳴数値シミュレーション装置。
前記数値計算部は、
前記磁化の更新後の数値を算出するために、前記磁化の更新前の数値を参照して前記磁化更新式を演算し、
前記算出対象方向での前記空間方向の前記偏微分の更新後の数値を算出するために、前記磁化の更新前の数値及び前記算出対象方向での前記空間方向の前記偏微分の更新前の数値を参照して前記空間偏微分更新式を演算し、
前記角周波数方向の前記偏微分の更新後の数値を算出するために、前記磁化の更新前の数値及び前記角周波数方向の前記偏微分の更新前の数値を参照して前記角周波数偏微分更新式を演算する、
請求項２記載の磁気共鳴数値シミュレーション装置。
【請求項4】
前記数値計算部は、
前記磁化の更新前の数値を参照して前記磁化更新式のうちの回転項を演算して磁化回転項値を算出し、
前記磁化回転項値に前記磁化更新式のうちの緩和項を適用して前記磁化の更新後の数値を算出する、
請求項３記載の磁気共鳴数値シミュレーション装置。
前記磁化の更新前の数値を参照して前記磁化更新式のうちの回転項を演算して磁化回転項値を算出し、
前記磁化回転項値に前記磁化更新式のうちの緩和項を適用して前記磁化の更新後の数値を算出する、
請求項３記載の磁気共鳴数値シミュレーション装置。
【請求項5】
前記数値計算部は、
前記磁化の更新前の数値及び前記算出対象方向での前記偏微分の更新前の数値を参照して前記空間偏微分更新式及び／又は前記角周波数偏微分更新式のうちの回転項を演算して偏微分回転項値を算出し、
前記偏微分回転項値に前記空間偏微分更新式及び／又は前記角周波数偏微分更新式のうちの緩和項を適用して前記算出対象方向での前記偏微分の更新後の数値を算出する、
請求項４記載の磁気共鳴数値シミュレーション装置。
前記磁化の更新前の数値及び前記算出対象方向での前記偏微分の更新前の数値を参照して前記空間偏微分更新式及び／又は前記角周波数偏微分更新式のうちの回転項を演算して偏微分回転項値を算出し、
前記偏微分回転項値に前記空間偏微分更新式及び／又は前記角周波数偏微分更新式のうちの緩和項を適用して前記算出対象方向での前記偏微分の更新後の数値を算出する、
請求項４記載の磁気共鳴数値シミュレーション装置。
【請求項6】
前記空間偏微分更新式は、互いに分離された、スライス選択傾斜磁場方向に対する直交方向の有効磁場ベクトル成分に依存する第１項と、スライス選択傾斜磁場方向の有効磁場ベクトル成分に依存する第２項と、前記有効磁場ベクトル成分に依存しない第３項とを有する、請求項５記載の磁気共鳴数値シミュレーション装置。
【請求項7】
前記数値計算部は、
ＲＦパルスの印加期間において、前記第１項、前記第２項及び前記第３項を演算し、
ＲＦパルスの非印加期間において、前記第２項及び前記第３項を演算する、
請求項６記載の磁気共鳴数値シミュレーション装置。
ＲＦパルスの印加期間において、前記第１項、前記第２項及び前記第３項を演算し、
ＲＦパルスの非印加期間において、前記第２項及び前記第３項を演算する、
請求項６記載の磁気共鳴数値シミュレーション装置。
【請求項8】
前記数値計算部は、ＲＦパルスの印加期間及び非印加期間において、前記第１項、前記第２項及び前記第３項を演算する、請求項６記載の磁気共鳴数値シミュレーション装置。
【請求項9】
前記数値計算部は、
前記算出対象方向での空間偏微分の更新後の数値を算出するために、前記磁化の更新前の数値及び前記算出対象方向での前記空間偏微分の更新前の数値を参照して前記第１項のうちの回転項を演算して第１回転項値を算出し、前記磁化の更新前の数値及び前記算出対象方向での前記空間偏微分の更新前の数値を参照して前記第２項のうちの回転項を演算して第２回転項値を算出し、前記磁化の更新前の数値及び前記算出対象方向での前記空間偏微分の更新前の数値を参照して前記第３項のうちの回転項を演算して第３回転項値を算出し、
前記第１回転項値と前記第２回転項値と前記第３回転項値との加算値に、前記空間偏微分更新式のうちの緩和項を適用して前記算出対象方向での前記空間偏微分の更新後の数値を算出する、
請求項６記載の磁気共鳴数値シミュレーション装置。
前記算出対象方向での空間偏微分の更新後の数値を算出するために、前記磁化の更新前の数値及び前記算出対象方向での前記空間偏微分の更新前の数値を参照して前記第１項のうちの回転項を演算して第１回転項値を算出し、前記磁化の更新前の数値及び前記算出対象方向での前記空間偏微分の更新前の数値を参照して前記第２項のうちの回転項を演算して第２回転項値を算出し、前記磁化の更新前の数値及び前記算出対象方向での前記空間偏微分の更新前の数値を参照して前記第３項のうちの回転項を演算して第３回転項値を算出し、
前記第１回転項値と前記第２回転項値と前記第３回転項値との加算値に、前記空間偏微分更新式のうちの緩和項を適用して前記算出対象方向での前記空間偏微分の更新後の数値を算出する、
請求項６記載の磁気共鳴数値シミュレーション装置。
【請求項10】
前記更新式は、複数の時間ステップ毎の前記磁化の時間発展と前記偏微分の時間発展とを記述する遷移行列を複数個組み合わせた結合遷移行列である、請求項２記載の磁気共鳴数値シミュレーション装置。
【請求項11】
前記結合遷移行列のうちの磁気共鳴の時間発展の演算に寄与する第１行列要素は、非ゼロ値を有し、
前記結合遷移行列のうちの磁気共鳴の時間発展の演算に寄与しない第２行列要素は、ゼロ値を有し、
前記数値計算部は、前記ゼロ値を有する前記第２行列要素の演算を実行しない、
請求項１０記載の磁気共鳴数値シミュレーション装置。
前記結合遷移行列のうちの磁気共鳴の時間発展の演算に寄与しない第２行列要素は、ゼロ値を有し、
前記数値計算部は、前記ゼロ値を有する前記第２行列要素の演算を実行しない、
請求項１０記載の磁気共鳴数値シミュレーション装置。
【請求項12】
前記結合遷移行列を記憶する記憶部を更に備え、
前記記憶部は、前記第１行列要素を記憶し、前記第２行列要素を記憶しない、
請求項１１記載の磁気共鳴数値シミュレーション装置。
前記記憶部は、前記第１行列要素を記憶し、前記第２行列要素を記憶しない、
請求項１１記載の磁気共鳴数値シミュレーション装置。
【請求項13】
既定の磁化の数値と前記空間方向各々での磁化の空間偏微分の数値とを含む入力ベクトルを前記更新式に適用して得られた出力ベクトルに基づいて前記結合遷移行列の各行列要素を算出する行列算出部を更に備える、請求項１０記載の磁気共鳴数値シミュレーション装置。
【請求項14】
パルスシーケンスのシーケンス断片毎に結合遷移行列を記憶する記憶部と、
演算対象のシーケンス断片に対応する結合遷移行列が前記記憶部に記憶されているか否かを判定する判定部と、を更に備え、
前記数値計算部は、
前記判定部により記憶されていると判定された場合、前記演算対象のシーケンス断片に対応する結合遷移行列を前記記憶部から読み出して使用し、
前記判定部により記憶されていないと判定された場合、前記演算対象のシーケンス断片に対応する結合遷移行列を算出し、前記算出された結合遷移行列を使用する、
請求項１０記載の磁気共鳴数値シミュレーション装置。
演算対象のシーケンス断片に対応する結合遷移行列が前記記憶部に記憶されているか否かを判定する判定部と、を更に備え、
前記数値計算部は、
前記判定部により記憶されていると判定された場合、前記演算対象のシーケンス断片に対応する結合遷移行列を前記記憶部から読み出して使用し、
前記判定部により記憶されていないと判定された場合、前記演算対象のシーケンス断片に対応する結合遷移行列を算出し、前記算出された結合遷移行列を使用する、
請求項１０記載の磁気共鳴数値シミュレーション装置。
【請求項15】
前記数値計算部は、Ａ／Ｄ変換器の予測出力値を算出する場合、スライス選択傾斜磁場方向の磁化成分に対する前記算出対象方向での前記偏微分の更新後の数値を算出せず、前記スライス選択傾斜磁場方向以外の傾斜磁場方向の磁化成分に対する前記算出対象方向での前記偏微分の更新後の数値を算出する、請求項１記載の磁気共鳴数値シミュレーション装置。
【請求項16】
アイソクロマットについて、磁化の更新前の数値と前記磁化に対する空間方向及び／又は角周波数方向各々での偏微分の更新前の数値とを入力し、
前記磁化の更新前の数値と前記偏微分の更新前の数値とのうちの一部の数値を使用して、磁気共鳴の時間発展を表す更新式を演算し、前記アイソクロマットについて、前記磁化の更新後の数値と前記磁化に対する前記空間方向及び／又は前記角周波数方向のうちの算出対象方向での前記偏微分の更新後の数値とを算出する、
こと具備する磁気共鳴数値シミュレーション方法。
前記磁化の更新前の数値と前記偏微分の更新前の数値とのうちの一部の数値を使用して、磁気共鳴の時間発展を表す更新式を演算し、前記アイソクロマットについて、前記磁化の更新後の数値と前記磁化に対する前記空間方向及び／又は前記角周波数方向のうちの算出対象方向での前記偏微分の更新後の数値とを算出する、
こと具備する磁気共鳴数値シミュレーション方法。
【請求項17】
コンピュータに、
アイソクロマットについて、磁化の更新前の数値と前記磁化に対する空間方向及び／又は角周波数方向各々での偏微分の更新前の数値とを入力させる機能と、
前記磁化の更新前の数値と前記偏微分の更新前の数値とのうちの一部の数値を使用して、磁気共鳴の時間発展を表す更新式を演算し、前記アイソクロマットについて、前記磁化の更新後の数値と前記磁化に対する前記空間方向及び／又は前記角周波数方向のうちの算出対象方向での前記偏微分の更新後の数値とを算出させる機能と、
を実現させる磁気共鳴数値シミュレーションプログラム。
アイソクロマットについて、磁化の更新前の数値と前記磁化に対する空間方向及び／又は角周波数方向各々での偏微分の更新前の数値とを入力させる機能と、
前記磁化の更新前の数値と前記偏微分の更新前の数値とのうちの一部の数値を使用して、磁気共鳴の時間発展を表す更新式を演算し、前記アイソクロマットについて、前記磁化の更新後の数値と前記磁化に対する前記空間方向及び／又は前記角周波数方向のうちの算出対象方向での前記偏微分の更新後の数値とを算出させる機能と、
を実現させる磁気共鳴数値シミュレーションプログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本明細書及び図面に開示の実施形態は、磁気共鳴数値シミュレーション装置、方法及びプログラムに関する。
【背景技術】
【0002】
磁気共鳴数値シミュレーションでは、ブロッホ方程式を用いて磁気共鳴数値イメージングをシミュレートする。ブロッホ方程式は磁場中の磁化の時間発展を記述する方程式であり、磁気共鳴イメージングの基礎方程式である。ブロッホ方程式の数値的解法としてＭＡＧＳＩ（Magnetic Gradient Simulation of Intra-voxel dephasing)がある。ＭＡＧＳＩは、ボクセルの中心に仮想分子を設定し、ボクセル内の磁場が線形であるとの仮定のもと、磁化の空間偏微分を用いてボクセル内の任意の位置の磁化を計算する。ただしＭＡＧＳＩは、磁化の空間偏微分を愚直に計算しているため、計算効率に改善の余地がある。
【先行技術文献】
【特許文献】
【0003】
【特許文献1】特開２０２３－００３３８９号公報
【非特許文献】
【0004】
【非特許文献1】Hidenori Takeshima. “A Fast and Practical Computation Method for Magnetic Resonance Simulators” Magnetic Resonance in Medicine, DOI： 10.1002/mrm.29646
【非特許文献2】Thies H. Jochimsen, Andreas Schafer, Roland Bammer, Michael E. Moseley, “Efficient simulation of magnetic resonance imaging with Bloch-Torrey equations using intra-voxel magnetization gradients” Journal of Magnetic Resonance, 2006; 180:29-38.
【非特許文献3】Graf C, Rund A, Aigner CS, Stollberger R. “Accuracy and Performance Analysis for Bloch and Bloch-McConnell Simulation Methods” Journal of Magnetic Resonance. 2021; 329: 107011.
【非特許文献4】Latta P, Gruwel MLH, Jellus V, Tomanek B. “Bloch simulations with intra-voxel spin dephasing” Journal of Magnetic Resonance. 2010; 203:44-51.
【発明の概要】
【発明が解決しようとする課題】
【0005】
本明細書及び図面に開示の実施形態が解決しようとする課題の一つは、磁気共鳴数値シミュレーションの計算効率を向上することである。ただし、本明細書及び図面に開示の実施形態により解決しようとする課題は上記課題に限られない。後述する実施形態に示す各構成による各効果に対応する課題を他の課題として位置づけることもできる。
【課題を解決するための手段】
【0006】
実施形態に係る磁気共鳴数値シミュレーション装置は、入力部と数値計算部とを有する。前記入力部は、アイソクロマットについて、磁化の更新前の数値と前記磁化に対する空間方向及び／又は角周波数方向各々での偏微分の更新前の数値とを入力する。前記数値計算部は、前記磁化の更新前の数値と前記偏微分の更新前の数値とのうちの一部の数値を使用して、磁気共鳴の時間発展を表す更新式を演算し、前記アイソクロマットについて、前記磁化の更新後の数値と前記磁化に対する前記空間方向及び／又は前記角周波数方向のうちの算出対象方向での前記偏微分の更新後の数値とを算出する。
【図面の簡単な説明】
【0007】
【発明を実施するための形態】
【0008】
以下、図面を参照しながら、本実施形態に係る磁気共鳴数値シミュレーション装置、方法及びプログラムについて詳細に説明する。
【0009】
本実施形態に係る磁気共鳴数値シミュレーション装置は、磁気共鳴数値シミュレーションを実行するコンピュータである。本実施形態に係る磁気共鳴数値シミュレーションは、検査対象モデルに設定された各ボクセル内のアイソクロマットに対する、所与のパルスシーケンスによる磁気共鳴の時間発展をシミュレーションする。具体的には、ボクセル内の磁場が線形であるとの仮定のもと、磁化の空間偏微分を用いてボクセル内の任意の位置の磁化を計算するＭＡＧＳＩが使用される。本実施形態に係る磁気共鳴数値シミュレーション装置、方法及びプログラムは、磁化及び／又は偏微分の計算効率を高めて計算効率を高めることを目的としている。
【0010】
（第１実施形態）
図１は、第１実施形態に係る磁気共鳴数値シミュレーション装置１の構成例を示す図である。図１に示すように、磁気共鳴数値シミュレーション装置１は、処理回路１１、通信機器１２、表示機器１３、入力機器１４及び記憶装置１５を有するコンピュータである。処理回路１１、通信機器１２、表示機器１３、入力機器１４及び記憶装置１５間のデータ通信は、バスを介して行われる。
図１は、第１実施形態に係る磁気共鳴数値シミュレーション装置１の構成例を示す図である。図１に示すように、磁気共鳴数値シミュレーション装置１は、処理回路１１、通信機器１２、表示機器１３、入力機器１４及び記憶装置１５を有するコンピュータである。処理回路１１、通信機器１２、表示機器１３、入力機器１４及び記憶装置１５間のデータ通信は、バスを介して行われる。
【0011】
処理回路１１は、ＣＰＵ（Central Processing Unit）等のプロセッサを有する。当該プロセッサが記憶装置１５等にインストールされた磁気共鳴数値シミュレーションプログラムを起動することにより、入力機能１１１、数値計算機能１１２及び出力制御機能１１３等を実現する。各機能１１１～１１３は単一の処理回路で実現される場合に限らない。複数の独立したプロセッサを組み合わせて処理回路を構成し、各プロセッサが磁気共鳴数値シミュレーションプログラムを実行することにより各機能１１１～１１３を実現するものとしても構わない。また、各機能１１１～１１３は、磁気共鳴数値シミュレーションプログラムを構成するモジュールとして実装されてもよいし、個別のハードウェアとして実装されてもよい。
【0012】
入力機能１１１により、処理回路１１は、磁気共鳴数値シミュレーションに用いられる種々の情報を入力する。一例として、処理回路１１は、計算対象のパルスシーケンスのシーケンス断片を入力する。シーケンス断片各々について数値計算機能１１２により磁気共鳴数値計算が行われる。入力機能１１１により、処理回路１１は、アイソクロマット（isochoromat）について、磁化の更新前の数値と当該磁化に対する空間方向及び角周波数方向各々での偏微分の更新前の数値とを入力する。アイソクロマットは、ボクセル内に設定された、同一周波数を有する核磁気モーメントの集合を意味する。更新前の数値としては、入力機器１４により手動的に入力された数値、数値計算機能１１２により算出された更新後の数値、予め設定された任意値が入力されればよい。
【0013】
数値計算機能１１２により、処理回路１１は、入力機能１１１により入力されたシーケンス断片各々について磁化ベクトルの更新演算を実行する。一例として、処理回路１１は、入力機能１１１により入力された磁化の更新前の数値と偏微分の更新前の数値とのうちの一部の数値を使用して、磁気共鳴の時間発展を表す更新式（以下、単に更新式と呼ぶ）を演算し、アイソクロマットについて、磁化の更新後の数値と磁化に対する空間方向及び／又は角周波数方向のうちの算出対象方向での空間偏微分の更新後の数値とを算出する。具体的には、処理回路１１は、偏微分の更新前の数値を使用せず、磁化の更新前の数値を使用して、更新式を演算し、磁化の更新後の数値を算出する。また、処理回路１１は、空間方向及び／又は角周波数方向のうちの非算出対象方向での偏微分の更新前の数値を使用せず、算出対象方向での偏微分の更新前の数値を使用して、更新式を演算し、算出対象方向での偏微分の更新後の数値を算出する。以下、磁化ベクトルのうちの磁化の数値を磁化値、磁化の偏微分の数値を偏微分値と呼ぶ。更新式としては、ブロッホ方程式又はその発展形が使用される。本実施形態に係るブロッホ方程式は、オリジナルのブロッホ方程式だけでなく、Torreyによる拡散（diffusion）項やMcConnellによる磁化移動（magnetization transfer）項を追加した発展形も含む。
【0014】
数値計算機能１１２により、処理回路１１は、更新後の磁化値と算出対象方向での更新後の偏微分値とに基づいて、磁気共鳴数値シミュレーションの目的に対応するデータ（以下、シミュレーションデータ）を生成してもよい。シミュレーションデータとしては、Ａ／Ｄ変換器の予測出力値や当該予測出力値に基づくＭＲ画像等が挙げられる。なお、磁化値と空間偏微分値とがシミュレーションデータに設定されてもよい。
【0015】
出力制御機能１１３により、処理回路１１は、数値計算機能１１２により算出された更新後の磁化値や偏微分値等のシミュレーションデータを出力する。一例として、処理回路１１は、シミュレーションデータを表示機器１３に表示したり、記憶装置１５に記憶したり、通信機器１２を介してほかのコンピュータに送信する。
【0016】
通信機器１２は、ＬＡＮ（Local Area Network）等を介して磁気共鳴イメージング装置、ワークステーション、ＰＡＣＳ（Picture Archiving and Communication System）、ＨＩＳ（Hospital Information System）、ＲＩＳ（Radiology Information System）等に接続するインタフェースである。通信機器１２は、各種データを接続先との間で送受信する。
【0017】
表示機器１３は、処理回路１１の出力制御機能１１３に従い各種データを表示する。表示機器１３としては、液晶ディスプレイ（ＬＣＤ：Liquid Crystal Display）、ＣＲＴ（Cathode Ray Tube）ディスプレイ、有機ＥＬディスプレイ（ＯＥＬＤ：Organic Electro Luminescence Display）、プラズマディスプレイ又は他の任意のディスプレイが適宜使用可能である。また、表示機器１３は、プロジェクタでもよい。
【0018】
入力機器１４は、ユーザからの各種の入力操作を受け付け、受け付けた入力操作を電気信号に変換して処理回路１１に出力する。具体的には、入力機器１４は、マウス、キーボード、トラックボール、スイッチ、ボタン、ジョイスティック、タッチパッド及びタッチパネルディスプレイ等の入力機器に接続されている。入力機器１４は、当該入力機器１４への入力操作に応じた電気信号を処理回路１１へ出力する。また、入力機器１４は、ネットワーク等を介して接続された他のコンピュータに設けられた入力機器でもよい。入力機器１４は、マイクロフォンにより収集された音声信号を指示信号に変換する音声認識装置でもよい。
【0019】
記憶装置１５は、各種データを記憶するＲＯＭ（Read Only Memory）やＲＡＭ（Random Access Memory）、ＨＤＤ（Hard Disk Drive）、ＳＳＤ（Solid State Drive）、集積回路記憶装置等の記憶装置である。記憶装置１５は、例えば、磁気共鳴数値シミュレーションプログラム等を記憶する。記憶装置１５は、上記記憶装置以外にも、ＣＤ（Compact Disc）、ＤＶＤ（Digital Versatile Disc）、フラッシュメモリ等の可搬型記憶媒体や、半導体メモリ素子等との間で種々の情報を読み書きする駆動装置であってもよい。記憶装置１５は、磁気共鳴数値シミュレーション装置１にネットワークを介して接続された他のコンピュータ内にあってもよい。
【0020】
上記の説明の通り、本実施形態に係る磁気共鳴数値シミュレーション１装置は、処理回路１１を有する。処理回路１１は、アイソクロマットについて、磁化の更新前の数値と磁化に対する空間方向及び／又は角周波数方向各々での偏微分の更新前の数値とを入力する。処理回路１１は、磁化の更新前の数値と偏微分の更新前の数値とのうちの一部の数値を使用して、磁気共鳴の時間発展を表す更新式を演算し、アイソクロマットについて、磁化の更新後の数値と空間方向及び／又は角周波数方向のうちの算出対象方向での偏微分の更新後の数値とを算出する。
【0021】
具体的には、処理回路１１は、更新前の磁化値を算出する場合、一部の数値として更新前の磁化値を使用する。換言すれば、処理回路１１は、一部の数値として更新前の偏微分値を使用しなくてもよい。更新前の算出対象方向での偏微分値を算出する場合、処理回路１１は、一部の数値として更新前の磁化値と当該算出対象方向での偏微分値を使用する。換言すれば、処理回路１１は、一部の数値として更新前の非算出対象方向での偏微分値を使用しなくてもよい。なお、空間方向での偏微分値及び／又は角周波数方向での偏微分値各々は、１個又は複数個の成分が存在し得る。例えば、空間方向での偏微分値は、一例として、リードアウト傾斜磁場方向での偏微分値、位相エンコード傾斜磁場方向での偏微分値、スライス選択傾斜磁場方向での偏微分値等の複数個の成分が存在し得る。算出対象方向での偏微分値は、これら複数個の成分のうちの一部を指すものとする。角周波数方向の偏微分値も同様、１個又は複数個の成分が存在し得る。すなわち、更新前の算出対象方向での偏微分値を算出する場合、処理回路１１は、空間方向及び／又は角周波数方向のうちの算出対象方向での偏微分の更新後の数値を算出するために、更新前の磁化値、算出対象方向の偏微分値及び非算出対象方向の偏微分値の全ての数値を使用する必要はなく、更新前の磁化値及び算出対象方向の偏微分値を使用すればよい。
【0022】
上記の構成によれば、本実施形態は、更新後の磁化値と偏微分値との全ての数値を使用する場合に比して、更新式の演算量を低減することができるので、結果的に、計算時間やメモリ消費量を低減することができ、ひいては、磁気共鳴数値シミュレーションの計算効率を向上することができる。
【0023】
ここで、本実施形態を非特許文献４に対して比較する。非特許文献４は、ブロッホ方程式におけるｘ、ｙ、ｚの各方向に対し、１０－６程度の微小なずれを与えたときの磁化を求めておき、各軸に対する磁化差分を微分とみなして位相を求めることでintra-voxel spin dephasingをシミュレートする。非特許文献４に比して、本実施形態は、微分を差分に置き換えることはせず、空間偏微分値を算出している。そのうえで本実施形態は、異なる偏微分方向での空間偏微分同士に相互作用がないことを利用し、更新後の数値の算出に寄与する成分にのみ更新式を適用する。これにより本実施形態は、結果の正確性を疎かにすることなく短時間で更新後の数値を算出することが可能になる。
【0024】
以下、第１実施形態に係る磁気共鳴数値シミュレーション装置１の一例について説明する。なお、以下に示す磁化ベクトルや更新式等は一例であり、本実施形態はこれに限定されるものではない。
【0025】
図２は、本実施形態に係る磁気共鳴数値シミュレーションの計算対象である磁化ベクトルＭ（ｔ）を構成する成分を例示する図である。図２に示す磁化ベクトルＭ（ｔ）は、時刻ｔにおけるボクセル内のアイソクロマットの磁化モーメントのベクトル和を表す。具体的には、磁化ベクトルＭ（ｔ）は、Ｍｘ、Ｍｙ、Ｍｚ、∂ｘＭｘ、∂ｘＭｙ、∂ｘＭｚ、∂ｙＭｘ、∂ｙＭｙ、∂ｙＭｚ、∂ｚＭｘ、∂ｚＭｙ、∂ｚＭｚ、∂ωＭｘ、∂ωＭｙ、∂ωＭｚ及びＭ０等の成分に分解できる。Ｍｘ、Ｍｙ及びＭｚは、磁化の空間方向の成分を表す。具体的には、Ｍｘ、Ｍｙ及びＭｚは、それぞれｘ方向、ｙ方向及びｚ方向の磁化成分を表す。∂ｘＭｘ、∂ｘＭｙ、∂ｘＭｚ、∂ｙＭｘ、∂ｙＭｙ、∂ｙＭｚ、∂ｚＭｘ、∂ｚＭｙ及び∂ｚＭｚは、磁化に対する複数の空間方向各々での空間偏微分を表す。具体的には、∂ｘＭｘ、∂ｘＭｙ及び∂ｘＭｚは、それぞれ磁化成分Ｍｘ、Ｍｙ及びＭｚに対するｘ方向での偏微分成分を表す。∂ｙＭｘ、∂ｙＭｙ、∂ｙＭｚは、それぞれ磁化成分Ｍｘ、Ｍｙ及びＭｚに対するｙ方向での偏微分成分を表す。∂ｚＭｘ、∂ｚＭｙ及び∂ｚＭｚは、それぞれ磁化成分Ｍｘ、Ｍｙ及びＭｚに対するｚ方向での偏微分成分を表す。∂ωＭｘ、∂ωＭｙ、∂ωＭｚは、それぞれ磁化成分Ｍｘ、Ｍｙ及びＭｚに対する角周波数方向での偏微分成分を表す。Ｍ０は、熱平衡状態におけるアイソクロマットの磁化を表す。Ｔは転置を意味する。
【0026】
なお、本実施形態においてｘ軸、ｙ軸及びｚ軸は、任意の方向に設定可能であるが、計算の便宜のため、ｘ軸にリードアウト傾斜磁場方向を、ｙ軸に位相エンコード傾斜磁場方向を、ｚ軸にスライス選択傾斜磁場方向を設定するものとする。また、角周波数方向での偏微分成分は、考慮しても良いし、しなくても良い。
【0027】
処理回路１１は、検査対象モデルに設定された複数のボクセル各々について、更新前の磁化ベクトルＭ（ｔ）に更新式を適用して、所与の時間ステップΔｔ後の磁化ベクトルＭ（ｔ＋Δｔ）を更新後の磁化ベクトルとして算出する。この際、処理回路１１は、各ボクセル内の磁場が線形であるとの仮定のもと、磁化の空間偏微分を用いてボクセル内の任意の位置の磁化を計算する。時間微分と空間偏微分とにより、ボクセル内の任意の位置における更新後の磁化を算出することができる。
【0028】
図３は、比較例に係る磁化ベクトルＭ（ｔ＋Δｔ）の算出過程を模式的に示す図である。図３に示すように、比較例においては、磁化ベクトルＭ（ｔ）のための１個の更新式ＦＡＬＬを演算する。更新式ＦＡＬＬは、磁化ベクトルＭ（ｔ）に含まれる各成分Ｍｘ、Ｍｙ、Ｍｚ、∂ｘＭｘ、∂ｘＭｙ、∂ｘＭｚ、∂ｙＭｘ、∂ｙＭｙ、∂ｙＭｚ、∂ｚＭｘ、∂ｚＭｙ、∂ｚＭｚ、∂ωＭｘ、∂ωＭｙ、∂ωＭｚを導く１個の計算式を意味する。この場合、これら全成分Ｍｘ、Ｍｙ、Ｍｚ、∂ｘＭｘ、∂ｘＭｙ、∂ｘＭｚ、∂ｙＭｘ、∂ｙＭｙ、∂ｙＭｚ、∂ｚＭｘ、∂ｚＭｙ、∂ｚＭｚ、∂ωＭｘ、∂ωＭｙ、∂ωＭｚ及びＭ０を使用して、更新式ＦＡＬＬを演算し、更新後の磁化値Ｍｘ，ｙ，ｚ（ｔ＋Δｔ）と更新後の空間偏微分値∂ｘＭｘ，ｙ，ｚ（ｔ＋Δｔ）、∂ｙＭｘ，ｙ，ｚ（ｔ＋Δｔ）及び∂ｚＭｘ，ｙ，ｚ（ｔ＋Δｔ）とが算出される。なおＭｘ，ｙ，ｚ（ｔ＋Δｔ）は、Ｍｘ（ｔ＋Δｔ）、Ｍｙ（ｔ＋Δｔ）、Ｍｚ（ｔ＋Δｔ）を縮約した表記である。
【0029】
図４は、磁化ベクトルＭ（ｔ）に関する算出対象と寄与項とを示す図である。算出対象は、磁化ベクトルを構成する成分のうちの着目している成分を意味する。空間偏微分値に関して言えば、算出対象は、ｘ方向、ｙ方向及びｚ方向等の偏微分する空間方向（以下、偏微分方向）により決定付けられる。寄与項は、算出対象を計算するために使用する成分、換言すれば、更新式ＦＡＬＬのうちの算出対象の計算に使用される項において使用される成分を意味する。寄与項以外の項（以下、非寄与項）は、算出対象を計算するために使用しない成分、換言すれば、更新式ＦＡＬＬのうちの算出対象の計算に使用される項に含まれない成分を意味する。空間偏微分値に関して言えば、非寄与項は、算出対象の偏微分方向以外の偏微分方向を意味する。
【0030】
図４の（１）に示すように、算出対象が磁化値Ｍｘ、Ｍｙ及びＭｚである場合、その寄与項は、Ｍｘ、Ｍｙ、Ｍｚ、Ｍ０であり、非寄与項は、全空間方向及び角周波数方向での偏微分値∂ｘＭｘ、∂ｘＭｙ、∂ｘＭｚ、∂ｙＭｘ、∂ｙＭｙ、∂ｙＭｚ、∂ｚＭｘ、∂ｚＭｙ、∂ｚＭｚ、∂ωＭｘ、∂ωＭｙ、∂ωＭｚである。図４の（２）に示すように、算出対象がｘ方向での空間偏微分値∂ｘＭｘ、∂ｘＭｙ、∂ｘＭｚである場合、その寄与項は、Ｍｘ、Ｍｙ、Ｍｚ、∂ｘＭｘ、∂ｘＭｙ、∂ｘＭｚであり、非寄与項は、ｙ方向、ｚ方向及びω方向での偏微分値∂ｙＭｘ、∂ｙＭｙ、∂ｙＭｚ、∂ｚＭｘ、∂ｚＭｙ、∂ｚＭｚ、∂ωＭｘ、∂ωＭｙ、∂ωＭｚとＭ０とである。図４の（３）に示すように、算出対象がｙ方向での空間偏微分値∂ｙＭｘ、∂ｙＭｙ、∂ｙＭｚである場合、その寄与項は、Ｍｘ、Ｍｙ、Ｍｚ、∂ｙＭｘ、∂ｙＭｙ、∂ｙＭｚであり、非寄与項は、ｘ方向、ｚ方向及びω方向での偏微分値∂ｘＭｘ、∂ｘＭｙ、∂ｘＭｚ、∂ｚＭｘ、∂ｚＭｙ、∂ｚＭｚ、∂ωＭｘ、∂ωＭｙ、∂ωＭｚとＭ０とである。図４の（４）に示すように、算出対象がｚ方向での空間偏微分値∂ｚＭｘ、∂ｚＭｙ及び∂ｚＭｚである場合、その寄与項は、Ｍｘ、Ｍｙ、Ｍｚ、∂ｚＭｘ、∂ｚＭｙ、∂ｚＭｚであり、非寄与項は、ｘ方向、ｙ方向及びω方向での偏微分値∂ｙＭｘ、∂ｙＭｙ、∂ｙＭｚ、∂ｚＭｘ、∂ｚＭｙ、∂ｚＭｚ、∂ωＭｘ、∂ωＭｙ、∂ωＭｚとＭ０とである。また、図４の（５）に示すように、算出対象がω方向での偏微分値∂ωＭｘ、∂ωＭｙ及び∂ωＭｚである場合、その寄与項は、Ｍｘ、Ｍｙ、Ｍｚ、∂ωＭｘ、∂ωＭｙ、∂ωＭｚであり、非寄与項は、全空間方向での空間偏微分値∂ｘＭｘ、∂ｘＭｙ、∂ｘＭｚ、∂ｙＭｘ、∂ｙＭｙ、∂ｙＭｚ、∂ｚＭｘ、∂ｚＭｙ、∂ｚＭｚとＭ０とである。
【0031】
なお、例えばBloch方程式にTorreyのDiffusion項が含まれる場合は、寄与項として全空間方向での空間偏微分値∂ｘＭｘ、∂ｘＭｙ、∂ｘＭｚ、∂ｙＭｘ、∂ｙＭｙ、∂ｙＭｚ、∂ｚＭｘ、∂ｚＭｙ、∂ｚＭｚが（Bloch方程式では非寄与項であったとしても）すべての算出対象について加わる。
【0032】
図４の（１）～（５）に示すように、磁化ベクトルの磁化成分、ｘ方向での偏微分成分、ｙ方向での偏微分成分、ｚ方向での偏微分成分及びω方向での偏微分成分の各々には寄与項と非寄与項とが存在することが分かる。例えば、磁化値を計算するために、全ての方向での偏微分値を使用しなくて良い。他の例では、ある方向での空間偏微分成分値を計算するために、他の空間方向での空間偏微分値を使用しなくて良い。この事は、異なる空間方向での偏微分成分同士は相互作用がないと表現することもできる。そこで、磁化値、ｘ方向偏微分値、ｙ方向偏微分値、ｚ方向偏微分値及びω方向偏微分値を分けて計算する。
【0033】
図５は、第１実施形態に係る磁化ベクトルの更新式の演算過程を模式的に示す図である。図５に示すように、第１実施形態においては、磁化ベクトルの更新式は、磁化成分の更新式（以下、磁化更新式）ＦＭと空間偏微分成分の更新式（以下、空間偏微分更新式）Ｆ∂ｐと角周波数偏微分成分の更新式（以下、角周波数偏微分更新式又はω方向偏微分更新式）Ｆ∂ωとに分けられる。空間偏微分更新式Ｆ∂ｐは、更にｘ方向での空間偏微分成分の更新式（以下、ｘ方向偏微分更新式）Ｆ∂ｘ、ｙ方向での空間偏微分成分の更新式（以下、ｙ方向偏微分更新式）Ｆ∂ｙ、ｚ方向での空間偏微分成分の更新式（以下、ｚ方向偏微分更新式）Ｆ∂ｚに分けられる。
【0034】
図５に示すように、処理回路１１は、磁化値Ｍｘ、ｙ、ｚ（ｔ＋Δｔ）を算出する場合、寄与項Ｍｘ、Ｍｙ、Ｍｚ、Ｍ０を使用し、非寄与項∂ｘＭｘ、∂ｘＭｙ、∂ｘＭｚ、∂ｙＭｘ、∂ｙＭｙ、∂ｙＭｚ、∂ｚＭｘ、∂ｚＭｙ、∂ｚＭｚ、∂ωＭｘ、∂ωＭｙ、∂ωＭｚを使用せず、磁化更新式ＦＭを演算する。処理回路１１は、ｘ方向偏微分値∂ｘＭｘ、ｙ、ｚ（ｔ＋Δｔ）を算出する場合、寄与項Ｍｘ、Ｍｙ、Ｍｚ、∂ｘＭｘ、∂ｘＭｙ、∂ｘＭｚ、Ｍ０を使用し、非寄与項∂ｙＭｘ、∂ｙＭｙ、∂ｙＭｚ、∂ｚＭｘ、∂ｚＭｙ、∂ｚＭｚ、∂ωＭｘ、∂ωＭｙ、∂ωＭｚを使用せず、ｘ方向偏微分更新式Ｆ∂ｘを演算する。処理回路１１は、ｙ方向偏微分値∂ｙＭｘ、ｙ、ｚ（ｔ＋Δｔ）を算出する場合、寄与項Ｍｘ、Ｍｙ、Ｍｚ、∂ｙＭｘ、∂ｙＭｙ、∂ｙＭｚ、Ｍ０を使用し、非寄与項∂ｘＭｘ、∂ｘＭｙ、∂ｘＭｚ、∂ｚＭｘ、∂ｚＭｙ、∂ｚＭｚ、∂ωＭｘ、∂ωＭｙ、∂ωＭｚを使用せず、ｙ方向偏微分更新式Ｆ∂ｙを演算する。処理回路１１は、ｚ方向偏微分値∂ｚＭｘ、ｙ、ｚ（ｔ＋Δｔ）を算出する場合、寄与項Ｍｘ、Ｍｙ、Ｍｚ、∂ｚＭｘ、∂ｚＭｙ、∂ｚＭｚ、Ｍ０を使用し、非寄与項∂ｘＭｘ、∂ｘＭｙ、∂ｘＭｚ、∂ｙＭｘ、∂ｙＭｙ、∂ｙＭｚ、∂ωＭｘ、∂ωＭｙ、∂ωＭｚを使用せず、ｚ方向偏微分更新式Ｆ∂ｚを演算する。処理回路１１は、ω方向偏微分値∂ωＭｘ、ｙ、ｚ（ｔ＋Δｔ）を算出する場合、寄与項Ｍｘ、Ｍｙ、Ｍｚ、∂ωＭｘ、∂ωＭｙ、∂ωＭｚ、Ｍ０を使用し、非寄与項∂ｘＭｘ、∂ｘＭｙ、∂ｘＭｚ、∂ｙＭｘ、∂ｙＭｙ、∂ｙＭｚ、∂ｚＭｘ、∂ｚＭｙ、∂ｚＭｚを使用せず、ω方向偏微分更新式Ｆ∂ωを演算する。
【0035】
上記の通り、第１実施形態に係る更新式ＦＡＬＬは、互いに分離された磁化更新式ＦＭと空間偏微分更新式Ｆ∂と角周波数偏微分更新式Ｆ∂ωとを有する。処理回路１１は、更新後の磁化値を算出するために、更新前の磁化値のみを参照して、磁化更新式ＦＭを演算する。ここで、「参照する」とは、記憶装置１５に記憶されている更新前の数値にアクセスして読み出すことを意味する。また、処理回路１１は、算出対象方向での更新後の空間偏微分値を算出するために、更新前の磁化値及び算出対象方向での空間偏微分値のみを参照して、空間偏微分更新式Ｆ∂ｐを演算する。空間偏微分更新式Ｆ∂ｐは、互いに分離されたｘ方向偏微分更新式Ｆ∂ｘ、ｙ方向偏微分更新式Ｆ∂ｙ及びｚ方向偏微分更新式Ｆ∂ｚを有する。処理回路１１は、更新後のｘ方向偏微分値を算出するために、更新前の磁化値及びｘ方向偏微分値のみを参照して、ｘ方向偏微分更新式Ｆ∂ｘを演算する。また、処理回路１１は、更新後のｙ方向偏微分値を算出するために、更新前の磁化値及びｙ方向偏微分値のみを参照して、ｙ方向偏微分更新式Ｆ∂ｙを演算する。また、処理回路１１は、更新後のｚ方向偏微分値を算出するために、更新前の磁化値及びｚ方向偏微分値のみを参照して、ｚ方向偏微分更新式Ｆ∂ｚを演算する。処理回路１１は、更新後の角周波数方向偏微分値を算出するために、更新前の磁化値及び角周波数方向偏微分値のみを参照して、角周波数偏微分更新式Ｆ∂ωを演算する。
【0036】
このように互いに相互作用しない成分に係る更新式同士を分離することにより、算出対象成分に係る更新式を演算する際、非寄与項を参照する必要がなくなる。非寄与項を参照せず、寄与項を参照して、算出対象成分に係る更新式を演算することにより、図３に示すような磁化ベクトルの更新式ＦＡＬＬを算出する場合に比して、更新式の計算時間や計算負荷を低減することが可能になる。
【0037】
次に、第１実施形態に係る磁気共鳴数値シミュレーションの処理手順について、実施例１と実施例２とに分けて説明する。
【0038】
［実施例１］
磁気共鳴の時間発展を表すブロッホ方程式は下記（１）式で表される。なお本実施形態において拡散（diffusion）項は考慮しないものとする。
磁気共鳴の時間発展を表すブロッホ方程式は下記（１）式で表される。なお本実施形態において拡散（diffusion）項は考慮しないものとする。
【0039】
【数1】
【0040】
Ｍは磁化ベクトル、Ｂは有効磁場ベクトル、γは磁気回転比である。Ｔ１は縦緩和時間、Ｔ２は横緩和時間である。上記ブロッホ方程式の各ｘ、ｙ、ｚ成分を展開するとそれぞれ下記（２）、（３）及び（４）式で表される。（２）、（３）及び（４）式は、それぞれｘ方向、ｙ方向及びｚ方向の磁化更新式に対応する。
【0041】
【数2】
【0042】
上記の（２）、（３）及び（４）式のように、磁化値の計算式には、磁化に関する項とＭ０に関する項が含まれているが、空間偏微分に関する項が含まれていない。このように、磁化値を計算する際に空間偏微分値を用いる必要がないため、磁化更新式と空間偏微分更新式とを分離することができる。分離することにより、分離しない場合に比して、更新後の磁化値を計算する時間を低減することが可能である。
【0043】
次に、ブロッホ方程式の空間方向偏微分の式を求める。ｘ、ｙ、ｚ方向、ｔ方向の微分はいずれも微小時間内で連続とする。ただし、磁場を区分的定数とした結果ｔが不連続になった位置は無視する。また、1つのアイソクロマットに対応した微小ボクセルに対し、ケミカルシフトも含めて静磁場を一定と考える。ｐ∈｛ｘ、ｙ、ｚ｝としてＭｘ（ｔ）、Ｍｙ（ｔ）及びＭｚ（ｔ）の空間偏微分はそれぞれ下記（５）、（６）及び（７）式で表される。
【0044】
【数3】
【0045】
上記の（５）、（６）及び（７）式のように、算出対象方向ｐでの空間偏微分値の計算式には、当該方向ｐでの偏微分項が含まれているが、他の方向ｑ（≠ｐ）での偏微分項が含まれていない。例えば、（５）、（６）及び（７）式においてｐ＝ｘとした場合、（５）、（６）及び（７）式において、算出対象方向ｘに関する項のみが存在し、非算出対象方向ｙ及び／又はｚを含む項は存在しない。すなわち、算出対象方向での空間偏微分値を計算する際に非算出対象方向での空間偏微分値を用いる必要がないため、算出対象方向毎に空間偏微分更新式を分離することができる。そのため、分離しない場合に比して、更新後の空間偏微分値を計算する時間を低減することが可能である。
【0046】
図６は、実施例１に係る磁気共鳴数値シミュレーションの処理手順を例示する図である。なお図６の開始時点において計算対象のパルスシーケンスが既に設定されているものとする。本実施形態に使用するパルスシーケンスは、特に制限はなく、如何なる種類のパルスシーケンスにも適用可能である。以下、一例として、本実施形態に使用するパルスシーケンスは、スピンエコー法によるパルスシーケンスであるとする。なお、図６に示す実施例では、角周波数方向偏微分については考慮しないものとする。
【0047】
図７は、スピンエコー法によるパルスシーケンスの一例を示す図である。図７において、「ＲＦ」のタイムラインは、送信回路からの駆動信号に応じて駆動するＲＦ送信コイルにより印加されるＲＦパルスのシーケンスを示す。「Ｇｘ」のタイムラインは、傾斜磁場電源からの駆動信号に応じて駆動するＧｘコイルにより印加されるリードアウト傾斜磁場パルスのシーケンスを示す。「Ｇｙ」のタイムラインは、傾斜磁場電源からの駆動信号に応じて駆動するＧｙコイルにより印加される位相エンコード傾斜磁場パルスのシーケンスを示す。「Ｇｚ」のタイムラインは傾斜磁場電源からの駆動信号に応じて駆動するＧｚコイルにより印加されるスライス選択傾斜磁場パルスのシーケンスを示している。図７に示すように、パルスシーケンスは、ＲＦパルスが印加されない期間である非印加期間４１と印加される期間である印加期間４２とに区分できる。
【0048】
図７に示すように、スピンエコー法では、９０°パルスの印加と共に第１のスライス選択傾斜磁場パルスの印加が行われ、１８０°パルスの印加と共に第２のスライス選択傾斜磁場パルスの印加が行われ、所与の磁場勾配の位相エンコードパルスが印加される。９０°パルスの印加からエコー時間の経過後に信号値のピークが来るように、図示しないエコー信号が発生する。エコー信号は、受信コイルにより受信され、受信されたアナログのエコー信号は、Ａ／Ｄ変換器（ＡＤＣ：Analog digital converter）によりデジタルのｋ空間データに変換される。ｋ空間データは、Ａ／Ｄ変換器の予測出力値の一例である。エコー信号の受信時において、リードアウト傾斜磁場パルスが印加される。９０°パルスの印加からリードアウト傾斜磁場パルスの印加までの時間がＴＲ（time of repetition）である。所与の位相エンコードステップ数に亘り上記パルスシーケンスが繰り返される。
【0049】
図６に示すように、処理回路１１は、入力機能１１１の実現により、シーケンス断片を入力する（ステップＳ１）。シーケンス断片は、任意のルールで分割されたパルスシーケンスの一部分である。パルスシーケンスは、所定時間毎に分割されてもよいし、「ＲＦ」「Ｇｘ」「Ｇｙ」「Ｇｚ」等のパルス種毎に分割されてもよいし、これらの組合せ単位で分割されてもよい。
【0050】
ステップＳ１が行われると処理回路１１は、数値計算機能１１２の実現により、磁化更新式及び空間偏微分更新式のパラメータ（以下、更新式パラメータ）を決定する（ステップＳ２）。ステップＳ２において処理回路１１は、ステップＳ１において入力されたシーケンス断片を利用して更新式パラメータを決定する。更新式パラメータは、検査対象であるモデルに配置された全てのアイソクロマット毎に決定される。更新式パラメータとしては、更新式に含まれるＴ１値やＴ２値、Ｂ値等が挙げられる。
【0051】
ステップＳ２が行われると処理回路１１は、入力機能１１１の実現により、更新前の磁化値Ｍ（ｔ）、空間偏微分値∂ｐＭ（ｔ）、熱平衡状態での磁化値Ｍ０を入力する（ステップＳ３）。初回の更新演算では、更新前のＭ（ｔ）及び∂ｐＭ（ｔ）は、任意の初期値が入力される。２回目以降の更新演算では、Ｍ（ｔ）及び∂ｐＭ（ｔ）については、前回の更新演算で得られたＭ（ｔ＋Δｔ）及び∂ｐＭ（ｔ＋Δｔ）が入力される。磁化値Ｍ０については任意の値が入力されればよい。更新前のＭ（ｔ）及び∂ｐＭ（ｔ）と磁化値Ｍ０とは、ステップＳ２以前に記憶装置１５に記憶されているものとする。処理回路１１は、記憶装置１５から更新前のＭ（ｔ）及び∂ｐＭ（ｔ）と磁化値Ｍ０とを読み出して自身に入力することとなる。
【0052】
ステップＳ３が行われると処理回路１１は、数値計算機能１１２の実現により、磁化値Ｍ（ｔ）及びＭ０のみを参照して、磁化更新式（２）式、（３）式及び（４）式を演算し、更新後の磁化値Ｍ（ｔ＋Δｔ）を算出する（ステップＳ４）。この際、処理回路１１は、空間偏微分値∂ｐＭ（ｔ）を参照する必要はない。これにより磁化値Ｍ（ｔ＋Δｔ）の計算速度が向上する。算出された更新後の磁化値Ｍ（ｔ＋Δｔ）は記憶装置１５に記憶される。磁化更新式はルンゲクッタ（Runge-kutta）法等を用いて演算することが可能である。
【0053】
ステップＳ４が行われると処理回路１１は、数値計算機能１１２の実現により、磁化値Ｍ（ｔ）及び∂ｐＭ（ｔ）のみを参照して、空間偏微分更新式（５）式、（６）式及び（７）式を演算し、更新後の空間偏微分値∂ｐＭ（ｔ＋Δｔ）を算出する（ステップＳ５）。ステップＳ５において処理回路１１は、算出対象方向ｐをｘ方向、ｙ方向及びｚ方向で順番に切り替えつつ各方向ｐについて空間偏微分更新式（５）式、（６）式及び（７）式を演算して算出対象方向ｐの空間偏微分値∂ｐＭ（ｔ＋Δｔ）を算出する。この際、処理回路１１は、非算出対象方向ｑの空間偏微分値∂ｑＭ（ｔ）を参照する必要はない。これにより各空間偏微分値∂ｐＭ（ｔ＋Δｔ）の計算速度が向上する。算出された更新後の空間偏微分値∂ｐＭ（ｔ＋Δｔ）は記憶装置１５に記憶される。空間偏微分更新式はルンゲクッタ法等を用いて演算することが可能である。
【0054】
ステップＳ３～Ｓ５は、検査対象モデルの全てのボクセルについて実行される。また処理回路１１は、各ボクセルについて、ボクセル内の磁場が線形であるとの仮定のもと、空間偏微分値を用いてボクセル内の任意の位置の磁化値を計算する。
【0055】
ステップＳ５が行われると処理回路１１は、数値計算機能１１２の実現により、更新演算を終了するか否かを判定する（ステップＳ６）。所定の更新終了条件を満たす場合に更新演算を終了すると判定される。更新終了条件は、例えば、更新回数が所定値に到達したこと等に設定されるとよい。更新演算を終了しないと判定された場合（ステップＳ６：ＮＯ）、処理回路１１は、ステップＳ６において更新演算を終了すると判定されるまで、次のタイミングステップについてステップＳ３からＳ６までの処理を反復する。
【0056】
更新演算を終了すると判定された場合（ステップＳ６：ＹＥＳ）、処理回路１１は、出力機能１１３の実現により、最終の反復回数における磁化値Ｍ（ｔ＋Δｔ）及び空間偏微分値∂ｐＭ（ｔ＋Δｔ）を出力する（ステップＳ７）。
【0057】
ステップＳ７が行われると処理回路１１は、数値計算機能１１２の実現により、シミュレーション終了を終了するか否かを判定する（ステップＳ８）。所定のシミュレーション終了条件を満たす場合にシミュレーションを終了すると判定される。シミュレーション終了条件は、例えば、計算対象のパルスシーケンスに含まれる全てのシーケンス断片について処理が完了したこと等に設定されるとよい。シミュレーション終了条件を満たさないと判定された場合（ステップＳ８：ＮＯ）、処理回路１１は、ステップＳ１に進み、入力機能１１１の実現により、次のシーケンス断片を入力する。そして処理回路１１は、当該シーケンス断片についてステップＳ２～Ｓ８を実行し、シミュレーション終了条件を満たすと判定されるまで、ステップＳ１～Ｓ８を繰り返す。これにより所定時刻毎の、パルスシーケンスによる検査対象モデルに含まれるボクセル毎且つアイソクロマット毎の磁化ベクトルＭ（ｔ）が得られることとなる。
【0058】
そしてシミュレーション終了条件を満たすと判定された場合（ステップＳ８：ＹＥＳ）、処理回路１１は、数値計算機能１１２の実現により、シミュレーション目的に対応するシミュレーションデータを算出する（ステップＳ９）。シミュレーションデータとしては、磁化ベクトルＭ（ｔ）やＩＱ信号、ＭＲ画像に設定されればよい。ＩＱ信号は、ＡＤＣにより収集されるｋ空間データを意味する。ＩＱ信号は、複数のアイソクロマットに対応する複数の磁化ベクトルＭ（Ｔ）の和ΣＭのうちｘ成分ΣＭＸとｙ成分ΣＭＹとの組合せを意味する。ＭＲ画像は、ＩＱ信号にフーリエ変換等の再構成演算を施すことにより生成される。
【0059】
ステップＳ９が行われると処理回路１１は、出力機能１１３の実現により、ステップＳ９において出力されたシミュレーションデータを出力する（ステップＳ１０）。処理回路１１は、シミュレーションデータを通信機器１２や表示機器１３、記憶装置１５等の任意の出力先に出力可能である。
【0060】
以上により、実施例１に係る磁気共鳴数値シミュレーションが終了する。
【0061】
［実施例２］
次に実施例２に係る磁気共鳴数値シミュレーションについて説明する。
次に実施例２に係る磁気共鳴数値シミュレーションについて説明する。
【0062】
実施例２に係る磁化更新式と空間偏微分更新式と角周波数偏微分更新式とは、回転項と緩和項とを含むものとする。実施例２に係る処理回路１１は、数値計算機能１１２の実現により、更新前の磁化値を参照して、磁化更新式のうちの回転項を演算して磁化回転項値を算出し、当該磁化回転項値に磁化更新式のうちの緩和項を適用して更新後の磁化値を算出する。同様に、処理回路１１は、数値計算機能１１２の実現により、更新前の磁化値及び算出対象方向での偏微分値を参照して空間偏微分更新式及び／又は角周波数偏微分更新式のうちの回転項を演算して偏微分回転項値を算出し、当該偏微分回転項値に空間偏微分更新式のうちの緩和項を適用して更新後の算出対象方向での偏微分値を算出する。
【0063】
空間偏微分更新式は、互いに分離された、スライス選択傾斜磁場方向に対する直交方向の有効磁場ベクトル成分に依存する第１項と、スライス選択傾斜磁場方向の有効磁場ベクトル成分に依存する第２項と、前記有効磁場ベクトル成分に依存しない第３項とを有する。実施例２に係る処理回路１１は、数値計算機能１１２の実現により、算出対象方向での更新後の空間偏微分値を算出するために、更新前の磁化値及び算出対象方向での空間偏微分値を参照して第１項のうちの回転項を演算して第１回転項値を算出し、更新前の磁化値及び算出対象方向での空間偏微分値を参照して第２項のうちの回転項を演算して第２回転項値を算出し、更新前の磁化値及び算出対象方向での空間偏微分値を参照して第３項のうちの回転項を演算して第３回転項値を算出する。そして処理回路１１は、第１回転項値と第２回転項値と第３回転項値との加算値に、空間偏微分更新式のうちの緩和項を適用して算出対象方向での更新後の空間偏微分値を算出する。
【0064】
まず、実施例２に係る更新方法の１例について具体的に説明する。なお、以下で示すすべての微分方程式の解は、例えばニュートン法やルンゲクッタ法を用い、微小時間を単位として繰り返し更新する形でも算出可能であり、微分方程式に基づく更新方法を示したすべての式は解法の１例である。
【0065】
上記（５）、（６）及び（７）式に示すブロッホ方程式の空間偏微分に対して、微小時間について、変数分離を用いて近似的に更新する。
【0066】
Ｂｘ（ｔ）及びＢｙ（ｔ）に依存する項、Ｂｚ（ｔ）に依存する項、それ以外の項に分けて順に更新すると考えると、例えば、微小時間内でのＭｘ成分の空間偏微分を表す（５）式は下記（８）、（９）及び（１０）式に分離できる。なお、Ｍｙ成分、Ｍｚ成分に対する式も同様に生成できる。
【0067】
【数4】
【0068】
空間的なＲＦパルスのゆらぎ∂Ｂｘ（ｔ）／∂ｐ、∂Ｂｙ（ｔ）／∂ｐをゼロとみなせば、下記（１１）式で表される。この微分方程式は回転をあらわし、ロドリゲス回転公式を用いて更新できる。Ｍｙ成分及びＭｚ成分についても同様である。（８）式又は（１１）式に基づく空間偏微分更新式は、Ｂｘ（ｔ）及びＢｙ（ｔ）に依存する項であり、上記「スライス選択傾斜磁場方向に対する直交方向の有効磁場ベクトル成分に依存する第１項」の一例である。以下、当該空間微分更新式をＢｘｙ依存の空間偏微分更新式と呼ぶ。
【0069】
【数5】
【0070】
（１１）式は（２）式に対し緩和を無視したロドリゲス回転公式に基づく非特許文献１記載の更新式と同様の式により、∂Ｍｘ、∂Ｍy、∂Ｍzをまとめたベクトルに対して、磁場Ｂx、Ｂｙ、Ｂzに対応する３ｘ３の回転行列を求め、それを乗ずることにより更新できる。
【0071】
Ｍｘ成分及びＭｙ成分に関する（９）式及び（１０）式に基づく空間偏微分更新式は、微小時間については、（１２）式で表される。なお、（１２）式のＢｚ（ｔ）バーは、微小時間Δｔ内において変動するＢｚ（ｔ）を特定の固定値とみなす場合の値を意味する。Ｍｚ成分に関する（９）式及び（１０）式に基づく空間偏微分更新式は、微小時間については、（１３）式で表される。
【数6】
【0072】
ここでＭｘｙ＝Ｍｘ＋ｉＭｙとする。ただし、（８）式がない場合には、Ｂｚ（ｔ）がｘ、ｙ、ｚ、ωに対する一次式で与えられる場合には微小時間でなくてもＭｘｙに対する更新式が適用できる。また、Ｍｚに対する更新式は微小時間でなくても適用できる。
【0073】
（１２）式により表される空間偏微分更新式は、Ｂｚ（ｔ）に依存する項であり、上記「スライス選択傾斜磁場方向の有効磁場ベクトル成分に依存する第２項」の一例である。以下、当該空間微分更新式をＢｚ依存の空間偏微分更新式と呼ぶ。（１３）式により表される空間偏微分更新式は、Ｂｘ（ｔ）、Ｂｙ（ｔ）、Ｂｚ（ｔ）に依存しない項であり、上記「有効磁場ベクトル成分に依存しない第３項」の一例である。以下、当該空間微分更新式をＢ非依存の空間偏微分更新式と呼ぶ。
【0074】
磁化ベクトルの更新式を、ロドリゲス回転公式を用いて表現すると、下記（１４）式で表すことが可能である。（１４）式は、非特許文献２の（５）式に相当する。
【0075】
【数7】
【0076】
（１４）式の右辺における磁化ベクトルＭ（ｔ）は、磁気共鳴数値シミュレーションにおける検査対象となるモデルに配置された複数の水素原子のうち、時刻ｔにおいて位置ｒにある領域に含まれるアイソクロマットの磁化ベクトルを表す。（１４）式の右辺におけるＲＲＦは、時間間隔Δｔにおいてフリップ角αを生じさせるＲＦパルスの印加によるｚ軸及びｘ軸についての回転項を示している。すなわち、回転項ＲＲＦは、磁化ベクトルＭ（ｔ）に対する回転の影響を示している。（１４）式の右辺におけるＲｒｅｌａｘは、時間間隔Δｔにおける磁化ベクトルＭ（ｔ）に対する磁気緩和の影響を示す緩和項である。（１４）式の右辺におけるＲｏｔｚ（θ）は、角度θに関するｚ軸についての回転項を示している。ここで、θｇは、時間間隔Δｔに亘る傾斜磁場の積分値に関連する角度である。また、θｉは、時間間隔Δｔに亘る静磁場の不均一性に伴う角度を示している。
【0077】
（１４）式をｘ方向、ｙ方向及びｚ方向各々で空間偏微分することにより、（５）、（６）及び（７）式にそれぞれ対応するｘ方向、ｙ方向及びｚ方向各々での空間偏微分更新式（ロドリゲス回転公式版の空間偏微分更新式）を得ることが可能である。ロドリゲス回転公式版のｘ方向偏微分更新式、ｙ方向偏微分更新式及びｚ方向偏微分更新式各々に変数分離を適用することにより、（８）式又は（１１）式に基づくＢｘｙ依存の空間偏微分更新式のロドリゲス回転公式版、（１２）式により表されるＢｚ依存の空間偏微分更新式のロドリゲス回転公式版、（１３）式により表されるＢ非依存の空間偏微分更新式のロドリゲス回転公式版を得ることが可能である。
【0078】
図８は、第１実施形態の実施例２に係る磁気共鳴数値シミュレーションの処理手順を例示する図である。図８は、図６のステップＳ３からステップＳ６までの間の処理を示している。ステップＳＡ１からＳＡ６までの各処理は、数値計算機能１１２により行われる。なお、図８に示す実施例では、角周波数方向偏微分については考慮しないものとする。
【0079】
ステップＳ３が行われると処理回路１１は、Ｍ（ｔ）及びＭ０のみを参照して、磁化更新式の回転項Ｒｏｔｚ（θｇ）、Ｒｏｔｚ（θｇ）及びＲＲＦを演算し、磁化回転項値を算出する（ステップＳＡ１）。Ｍ（ｔ）及びＭ０のみを参照することにより、磁化更新式の回転項Ｒｏｔｚ（θｇ）、Ｒｏｔｚ（θｇ）及びＲＲＦの演算を効率良く行うことができる。
【0080】
ステップＳＡ１が行われると処理回路１１は、ステップＳＡ１において算出された磁化回転項値に磁化更新式の緩和項Ｒｒｅｌａｘを適用して更新後の磁化値Ｍ（ｔ＋Δｔ）を算出する（ステップＳＡ２）。更新後の磁化値Ｍ（ｔ＋Δｔ）は記憶装置１５に記憶される。
【0081】
ステップＳＡ２が行われると処理回路１１は、Ｍ（ｔ）及び∂ｐＭ（ｔ）のみを参照してＢｘｙ依存の空間偏微分更新式の回転項Ｒｏｔｚ（θｇ）、Ｒｏｔｚ（θｇ）及びＲＲＦを演算し、第１回転項値を算出する（ステップＳＡ３）。Ｍ（ｔ）及び∂ｐＭ（ｔ）のみを参照することにより、Ｂｘｙ依存の空間偏微分更新式の回転項Ｒｏｔｚ（θｇ）、Ｒｏｔｚ（θｇ）及びＲＲＦの演算を効率良く行うことができる。
【0082】
ステップＳＡ３が行われると処理回路１１は、Ｍ（ｔ）及び∂ｐＭ（ｔ）のみを参照してＢｚ依存の空間偏微分更新式の回転項Ｒｏｔｚ（θｇ）、Ｒｏｔｚ（θｇ）及びＲＲＦを演算し、第２回転項値を算出する（ステップＳＡ４）。Ｍ（ｔ）及び∂ｐＭ（ｔ）のみを参照することにより、Ｂｚ依存の空間偏微分更新式の回転項Ｒｏｔｚ（θｇ）、Ｒｏｔｚ（θｇ）及びＲＲＦの演算を効率良く行うことができる。なお、Ｂｚ（ｔ）バーが未算出である場合、上記方法に従いＢｚ（ｔ）バーを算出した後に回転項Ｒｏｔｚ（θｇ）、Ｒｏｔｚ（θｇ）及びＲＲＦの演算を行うと良い。
【0083】
ステップＳＡ４が行われると処理回路１１は、Ｍ（ｔ）及び∂ｐＭ（ｔ）のみを参照してＢ非依存の空間偏微分更新式の回転項Ｒｏｔｚ（θｇ）、Ｒｏｔｚ（θｇ）及びＲＲＦを演算し、第３回転項値を算出する（ステップＳＡ５）。Ｍ（ｔ）及び∂ｐＭ（ｔ）のみを参照することにより、Ｂ非依存の空間偏微分更新式の回転項Ｒｏｔｚ（θｇ）、Ｒｏｔｚ（θｇ）及びＲＲＦの演算を効率良く行うことができる。
【0084】
ステップＳＡ５が行われると処理回路１１は、ステップＳＡ３において算出された第１回転項値とステップＳＡ４において算出された第２回転項値とステップＳＡ４において算出された第３回転項値とを加算する（ステップＳＡ６）。ステップＳＡ６により第１回転項値と第２回転項値と第３回転項値との加算値が算出される。
【0085】
ステップＳＡ６が行われると処理回路１１は、ステップＳＡ６で得られた加算値に偏微分更新式の緩和項Ｒｒｅｌａｘを適用して更新後の空間偏微分値∂ｐＭ（ｔ＋Δｔ）を算出する（ステップＳＡ７）。更新後の空間偏微分値∂ｐＭ（ｔ＋Δｔ）は記憶装置１５に記憶される。
【0086】
ステップＳＡ６が行われると、図６に示すステップＳ６が実行される。
【0087】
上記の通り、実施例２によれば、ロドリゲス回転公式を使用して記述した更新式においても、磁化更新式と空間偏微分更新式と各周波数偏微分更新式とを分離し、非寄与項を使用せず寄与項を使用して効率良く磁化値と空間偏微分値とを計算することが可能になる。また、空間偏微分更新式においても算出対象方向毎に空間偏微分更新式を分離することできるので、非寄与項を使用せず寄与項を使用して効率良く各算出対象方向での空間偏微分値を計算することが可能になる。
【0088】
なお、ステップＳＡ３、ＳＡ４及びＳＡ５の処理順序はこれに限定されず、如何なる順番に行われてもよい。また、必ずしもＢｘｙ依存の空間偏微分更新式（第１項）、Ｂｚ依存の空間偏微分更新式（第２項）及びＢ非依存の空間偏微分更新式（第３項）の全てを演算する必要はなく、必要に応じて演算対象が適宜選択されればよい。例えば、処理回路１１は、ＲＦパルスの印加期間において、第１項、第２項及び第３項を演算し、ＲＦパルスの非印加期間において、第２項及び第３項を演算してもよい。他の例として、処理回路１１は、ＲＦパルスの印加期間及び非印加期間の双方において、第１項、第２項及び第３項を演算してもよい。
【0089】
また、得られたｘ、ｙ、ｚ方向の偏微分を用いてＡＤＣサンプリングをシミュレートする場合は、例えば非特許文献２記載の手法にしたがい、単位直方体内に分子が一様に置かれているモデルを用い、Ｍ（ｔ）、∂xＭ（ｔ）、∂yＭ（ｔ）を用いて角度（∂／∂ｐ）φを非特許文献２の（１６）式によって算出し、そのｓｉｎｃ関数の積分値を乗ずることで磁場の空間変動をサンプリング値に反映できる。また、ω方向の偏微分を用いてＴ２スター減衰をシミュレートする場合は、非特許文献２記載の手法（特定の周波数領域に対してのみ一様に分子が置かれている）を用いても良いし、例えば別の手法として、Ｔ２スター減衰が周波数に対するローレンツ分布にしたがった確率で分子が置かれるモデルを用い、ω方向の偏微分に対してｅｘｐ｛（－１／Ｔ２´）・（∂／∂ω）φ）｝を乗じても良い。Ｔ２´は磁場不均一に起因するＴ２スター減衰（Ｔ２
＊）とＴ２減衰の違いをあらわし、１／Ｔ２´＝（１／Ｔ２
＊）－（１／Ｔ２）で与えられる。
【0090】
（第２実施形態）
第２実施形態に係る処理回路１１は、磁気共鳴の時間発展を記述する更新式として結合遷移行列（combined transition matrix）を使用する。結合遷移行列は、複数の時間ステップ毎の磁化の時間発展と偏微分の時間発展とを記述する遷移行列Ｔ（ｔ＋Δｔｋ－１）をＫ（Ｋ≧２）個組み合わせた行列である。ｋは１以上の自然数である。Ｋ個の遷移行列Ｔ（ｔ）、Ｔ（ｔ＋Δｔ１）、・・・、Ｔ（ｔ＋ΔｔＫ－１）を結合した結合遷移行列Ｔｃｏｍｂｉｎｅｄは、下記（１５）式で表現できる。Δｔｋは、ｋｔから（ｋ＋１）ｔまでの時間間隔を表す。
第２実施形態に係る処理回路１１は、磁気共鳴の時間発展を記述する更新式として結合遷移行列（combined transition matrix）を使用する。結合遷移行列は、複数の時間ステップ毎の磁化の時間発展と偏微分の時間発展とを記述する遷移行列Ｔ（ｔ＋Δｔｋ－１）をＫ（Ｋ≧２）個組み合わせた行列である。ｋは１以上の自然数である。Ｋ個の遷移行列Ｔ（ｔ）、Ｔ（ｔ＋Δｔ１）、・・・、Ｔ（ｔ＋ΔｔＫ－１）を結合した結合遷移行列Ｔｃｏｍｂｉｎｅｄは、下記（１５）式で表現できる。Δｔｋは、ｋｔから（ｋ＋１）ｔまでの時間間隔を表す。
【0091】
【数8】
【0092】
結合遷移行列Ｔｃｏｍｂｉｎｅｄを使用した磁化ベクトルＭ（ｔ）の更新式は、下記（１６）で表現できる。
【0093】
【数9】
【0094】
図４等に示すように、結合遷移行列Ｔｃｏｍｂｉｎｅｄは、磁化ベクトルの成分数の２乗個の要素を含み、多くの要素を含む。結合遷移行列Ｔｃｏｍｂｉｎｅｄの個々の要素をＴａ，ｂで表すと、結合遷移行列Ｔｃｏｍｂｉｎｅｄは下記（１７）式で表現できる。なお、（１７）式では、磁化ベクトルの成分数を「ｎ」とおくことにする。
【0095】
【数10】
【0096】
結合遷移行列を使用する場合、１ボクセル毎に磁化ベクトルの成分数の２乗個の要素を有する遷移行列のＫ回の乗算を繰り返すこととなる。検査対象モデルが３次元の場合、ボクセルの個数は１００万を超えることも多く、計算に伴う計算量やメモリ使用量等の計算負荷が著しい。
【0097】
第２実施形態に係る磁気共鳴数値シミュレーション装置は、第１実施形態に係る、寄与項のみを参照する方法を結合遷移行列に応用して、磁気共鳴数値シミュレーションの高速化を図る。以下、第２実施形態に係る磁気共鳴数値シミュレーション装置について説明する。なお以下の説明において、第１実施形態と略同一の機能を有する構成要素については、同一符号を付し、必要な場合にのみ重複説明する。
【0098】
図９は、第２実施形態に係る磁気共鳴数値シミュレーション装置１の構成例を示す図である。図９に示すように、磁気共鳴数値シミュレーション装置１の処理回路１１は、入力機能１１１、数値計算機能１１２及び出力制御機能１１３の他、判定機能１１４及び行列算出機能１１５を実現する。各機能１１１～１１５は単一の処理回路で実現される場合に限らない。複数の独立したプロセッサを組み合わせて処理回路を構成し、各プロセッサが物理シミュレーションプログラムを実行することにより各機能１１１～１１５を実現するものとしても構わない。また、各機能１１１～１１５は、磁気共鳴数値シミュレーションプログラムを構成するモジュールとして実装されてもよいし、個別のハードウェアとして実装されてもよい。
【0099】
数値計算機能１１２の実現により、処理回路１１は、更新前の磁化値及び空間偏微分値を結合遷移行列に適用して更新後の磁化値及び空間偏微分値を算出する。結合遷移行列は、予め行列算出機能１１５により算出され、記憶装置１５に記憶されているものとする。結合遷移行列の各行列要素の更新式パラメータがシーケンス断片毎に異なるので、結合遷移行列は、シーケンス断片毎に算出され記憶される。一例として、結合遷移行列は、シーケンス断片を表す識別子に関連付けて記憶されるとよい。
【0100】
判定機能１１４の実現により、処理回路１１は、演算対象のシーケンス断片に対応する結合遷移行列が記憶装置１５に記憶されているか否かを判定する。処理回路１１は、記憶されていると判定された場合、演算対象のシーケンス断片に対応する結合遷移行列を記憶装置１５から読み出して使用する。
【0101】
行列算出機能１１５の実現により、処理回路１１は、既定の磁化値と複数の空間方向各々での空間偏微分値とを含む入力ベクトルを更新式に適用して得られた出力ベクトルに基づいて結合遷移行列の各行列要素を算出する。結合遷移行列は、記憶装置１５に記憶される。一例として、処理回路１１は、判定機能１１４により演算対象のシーケンス断片に対応する結合遷移行列が記憶装置１５に記憶されていないと判定された場合、当該シーケンス断片に対応する結合遷移行列を算出して使用してもよい。なお、予めシーケンスがわかっている場合は、シミュレーションより前に結合遷移行列を予め計算しても良い。また、結合遷移行列の算出対象とするシーケンス断片は、例えばＲＦ照射の有無を基準として切り出すことができるが、この基準に限らずどのような基準で切り出したシーケンスであっても良い。また、すべてのシーケンス断片に対して結合遷移行列の算出を行う必要はなく、一部のシーケンス断片についてのみ結合遷移行列を計算しても良い。
【0102】
以下、第２実施形態に磁気共鳴数値シミュレーション装置１について詳細に説明する。
【0103】
まず、結合遷移行列に含まれる行列要素について説明する。
【0104】
結合遷移行列のうちの磁気共鳴の時間発展の演算に寄与する第１行列要素は、非ゼロ値を有し、前記結合遷移行列のうちの磁気共鳴の時間発展の演算に寄与しない第２行列要素は、ゼロ値を有する。具体的には、第１行列要素は、磁化ベクトルのある成分と当該成分に相互作用する成分との組合せに対応する行列要素であり、第２行列要素は、磁化ベクトルのある成分と当該成分に相互作用しない成分との組合せに対応する行列要素である。
【0105】
より詳細には、Ｍｘ、Ｍｙ、Ｍｚは、Ｍｘ、Ｍｙ、Ｍｚ及びＭ０以外の成分に相互作用しない。したがってＭｘ、Ｍｙ、Ｍｚ各々について、非ゼロ値の行列要素は、Ｍｘ、Ｍｙ、Ｍｚ及びＭ０の４個である。∂ｐＭは、∂ｑＭ（ｐ≠ｑ）及びＭ０に相互作用しない。したがって∂ｐＭｘ、∂ｐＭｙ、∂ｐＭｚ各々について、非ゼロ値の行列要素は、Ｍｘ、Ｍｙ、Ｍｚ、∂ｐＭｘ、∂ｐＭｙ、∂ｐＭｚの６個である。また、Ｍ０は定数であるから、Ｍ０に関する行は全て除外できる。
【0106】
第２実施形態に係る処理回路１１は、ゼロ値を有する第２行列要素の演算を実行しない。すなわち、ゼロ値を有する第２行列要素は結合遷移行列の計算式に組み込まれず、当該計算式は非ゼロ値を有する第１行列要素により構成される。これにより、ゼロ値を有する多くの第２行列要素の演算を省略できるので、計算速度が向上する。
【0107】
結合遷移行列は記憶装置１５に記憶されるが、記憶装置１５は、非ゼロ値を有する第１行列要素を記憶し、ゼロ値を有する第２行列要素を記憶しない。これにより結合遷移行列を保持するメモリ容量を削減することができる。
【0108】
次に、行列算出機能１１５による結合遷移行列の算出について説明する。
【0109】
拡散項がない場合、各∂ｐＭは∂ｑＭ（ｐ≠ｑ）に対して独立である。結合遷移行列は、入力ベクトルに対応する出力ベクトルに基づいて算出できる。入力ベクトルは、既定の磁化値と複数の空間方向各々での空間偏微分値とを含む。出力ベクトルは、入力ベクトルに複合遷移行列を適用することにより得られる磁化ベクトルである。なお、Ｍ０は各シミュレーションにおいて不変である。
【0110】
入力ベクトルは、すべての独立な単位ベクトルを用意し、第１実施形態で説明した演算を実行すれば結合遷移行列が容易に得られる。あるいは、結合遷移行列はゼロの要素が多いため、それらの計算の一部を省略できる。例えばω方向を考慮しない場合、下記の７個のみを計算しても結合遷移行列を得ることができる。なお、ω方向を考慮する場合、Ｖ４～Ｖ６に∂ωに対応した項を追加すればよい。
【0111】
Ｖ１．Ｍｘ＝１、残り０
Ｖ２．Ｍｙ＝１、残り０
Ｖ３．Ｍｚ＝１、残り０
Ｖ４．∂ｘＭｘ＝∂yＭx＝∂zＭx＝１、残り０
Ｖ５．∂ｘＭｙ＝∂ｙＭｙ＝∂ｚＭｙ＝１、残り０
Ｖ６．∂ｘＭｚ＝∂ｙＭｚ＝∂ｚＭｚ＝１、残り０
Ｖ７．全て０
Ｖ２．Ｍｙ＝１、残り０
Ｖ３．Ｍｚ＝１、残り０
Ｖ４．∂ｘＭｘ＝∂yＭx＝∂zＭx＝１、残り０
Ｖ５．∂ｘＭｙ＝∂ｙＭｙ＝∂ｚＭｙ＝１、残り０
Ｖ６．∂ｘＭｚ＝∂ｙＭｚ＝∂ｚＭｚ＝１、残り０
Ｖ７．全て０
【0112】
Ｖ１～Ｖ７各々について、第１実施形態で説明した演算を実行する。Ｍ0の更新式に対応する結合遷移行列の要素は既知である。ゼロの要素を利用すれば、Ｖ１、Ｖ２、Ｖ３、Ｖ７からＭｘ、Ｍｙ、Ｍｚ、Ｍ0の更新式に対応する結合遷移行列の要素が得られる。同様に、ゼロの要素を利用し、各∂pの相互干渉がないことを利用すれば、Ｖ１、Ｖ２、Ｖ３、Ｖ４～Ｖ６の∂pに対応する項から∂pに対応する結合遷移行列の要素が得られる。
【0113】
ロドリゲス回転行列を用いる場合は、第１実施形態の実施例２の図８に示すステップＳＡ１、ＳＡ３及びＳＡ４各々の計算について中間計算の一部を共有することにより、結合遷移行列を算出することが可能である。あるいは、行列を陽に書き下して算出することが可能である。初期値を利用してソフトウェア実装する場合は、アイソクロマットごとに７個の初期ベクトルに対するシミュレーションを一括実行することにより、結合遷移行列の算出に要する時間を短縮することが可能である。遷移行列の行列積を利用して結合遷移行列を求める場合は、各部分行列で乗算を分けることにより、結合遷移行列の算出に要する時間を短縮することが可能である。いずれの場合も、ｓｉｎ、ｃｏｓ、ｅｘｐの値を可能な限り共有することにより、結合遷移行列の算出に要する時間を短縮することが可能である。
【0114】
［実施例１］
図１０は、第２実施形態の実施例１に係る磁気共鳴数値シミュレーションの処理手順を例示する図である。図１０は、図６のステップＳ３からステップＳ５までの間の処理を示している。ステップＳＢ１からＳＢ２までの各処理は、数値計算機能１１２により行われる。
図１０は、第２実施形態の実施例１に係る磁気共鳴数値シミュレーションの処理手順を例示する図である。図１０は、図６のステップＳ３からステップＳ５までの間の処理を示している。ステップＳＢ１からＳＢ２までの各処理は、数値計算機能１１２により行われる。
【0115】
図１０に示すように、ステップＳ３が行われると処理回路１１は、ステップＳ１において入力されたシーケンス断片に対応する結合遷移行列Ｔｃｏｍｂｉｎｅｄ（ｔ）を記憶装置１５から読み出す（ステップＳＢ１）。
【0116】
ステップＳＢ１が行われると処理回路１１は、Ｍ（ｔ）、∂ｐＭ（ｔ）及びＭ０を、ステップＳＢ１において読み出した結合遷移行列Ｔｃｏｍｂｉｎｅｄ（ｔ）に適用して、更新後のＭ（ｔ＋Δｔ）、∂ｐＭ（ｔ＋Δｔ）を算出する（ステップＳＢ２）。ステップＳＢ２において処理回路１１は、ゼロ値の第２行列要素の演算は実行せず、非ゼロ値の第１行列要素を実行する。これにより結合遷移行列の計算時間を短縮することが可能である。
【0117】
ステップＳＢ２が行われると、図６に示すステップＳ６が実行される。
【0118】
【0119】
図１１に示すように、ステップＳ３が行われると処理回路１１は、判定機能１１４の実現により、ステップＳ１において入力されたシーケンス断片に対応する結合遷移行列Ｔｃｏｍｂｉｎｅｄ（ｔ）が記憶装置１５に記憶されているか否かを判定する（ステップＳＣ１）。
【0120】
ステップＳＣ１において結合遷移行列Ｔｃｏｍｂｉｎｅｄ（ｔ）が記憶されていないと判定された場合（ステップＳＣ１：ＮＯ）、処理回路１１は、行列算出機能１１５の実現により、ステップＳ１において入力されたシーケンス断片に対応する結合遷移行列Ｔｃｏｍｂｉｎｅｄ（ｔ）を算出する（ステップＳＣ２）。
【0121】
ステップＳＣ２が行われると処理回路１１は、出力制御機能１１３の実現により、ステップＳＣ２において算出された結合遷移行列Ｔｃｏｍｂｉｎｅｄ（ｔ）を記憶装置１５にキャッシュする（ステップＳＣ３）。結合遷移行列Ｔｃｏｍｂｉｎｅｄ（ｔ）は、ステップＳ１において入力されたシーケンス断片の識別子に関連付けて記憶される。
【0122】
ステップＳＣ３が行われた場合又はステップＳＣ１において結合遷移行列Ｔｃｏｍｂｉｎｅｄ（ｔ）が記憶されていると判定された場合（ステップＳＣ１：ＹＥＳ）、処理回路１１は、ステップＳ１において入力されたシーケンス断片に対応する結合遷移行列Ｔｃｏｍｂｉｎｅｄ（ｔ）を記憶装置１５から読み出す（ステップＳＣ４）。
【0123】
ステップＳＣ４が行われると処理回路１１は、数値計算機能１１２の実現により、Ｍ（ｔ）、∂ｐＭ（ｔ）及びＭ０を、ステップＳＢ１において読み出した結合遷移行列Ｔｃｏｍｂｉｎｅｄ（ｔ）に適用して、更新後のＭ（ｔ＋Δｔ）、∂ｐＭ（ｔ＋Δｔ）を算出する（ステップＳＣ５）。
【0124】
ステップＳＣ５が行われると、図６に示すステップＳ６が実行される。
【0125】
実施例２によれば、未知のシーケンス断片が検出される毎に結合遷移行列Ｔｃｏｍｂｉｎｅｄ（ｔ）を計算し、その結果を動的にキャッシュしながら再利用することで、磁気共鳴数値シミュレーションの高速化が達成される。
【0126】
（変形例）
上記の磁気共鳴数値シミュレーション装置１による処理は一例であり、本実施形態の要旨を変更しない限り種々の追加、変更及び／又は削除が可能である。
上記の磁気共鳴数値シミュレーション装置１による処理は一例であり、本実施形態の要旨を変更しない限り種々の追加、変更及び／又は削除が可能である。
【0127】
［変形例１］
上記実施形態に係る処理回路１１は、磁化の更新前の数値と偏微分の更新前の数値とのうちの寄与項に属する全ての更新前の数値を使用して、算出対象に関する更新式を演算し、算出対象の更新後の数値とを算出するものとした。しかしながら、本実施形態はこれに限定されない。処理回路１１は、寄与項に属する全ての更新前の数値を使用する必要はなくてもよく、寄与項に属する全ての更新前の数値のうちの一部の数値のみを使用してもよい。一例として、変形例１に係る処理回路１１は、寄与項に属する更新前の数値のうちのゼロ又は略ゼロと見做せる数値を使用しなくてもよい。
上記実施形態に係る処理回路１１は、磁化の更新前の数値と偏微分の更新前の数値とのうちの寄与項に属する全ての更新前の数値を使用して、算出対象に関する更新式を演算し、算出対象の更新後の数値とを算出するものとした。しかしながら、本実施形態はこれに限定されない。処理回路１１は、寄与項に属する全ての更新前の数値を使用する必要はなくてもよく、寄与項に属する全ての更新前の数値のうちの一部の数値のみを使用してもよい。一例として、変形例１に係る処理回路１１は、寄与項に属する更新前の数値のうちのゼロ又は略ゼロと見做せる数値を使用しなくてもよい。
【0128】
［変形例２］
変形例２に係る処理回路１１は、Ａ／Ｄ変換器（ＡＤＣ）の予測出力値を算出する場合、スライス選択傾斜磁場方向（ｚ方向）の磁化成分に対する算出対象方向での更新後の偏微分値を算出せず、スライス選択傾斜磁場方向以外の傾斜磁場方向の磁化成分に対する算出対象方向での更新後の偏微分値を算出する。例えば、変形例２に係る処理回路１１は、Ａ／Ｄ変換器（ＡＤＣ）の予測出力値を算出する場合、∂ｐＭｚを算出せず、ＩＱ信号に対応する∂ｐＭｘ及び∂ｐＭｙのみを算出すればよい。∂ｐＭｚを算出する計算時間を削減し、Ａ／Ｄ変換器（ＡＤＣ）の予測出力値の算出時間を短縮することが可能になる。
変形例２に係る処理回路１１は、Ａ／Ｄ変換器（ＡＤＣ）の予測出力値を算出する場合、スライス選択傾斜磁場方向（ｚ方向）の磁化成分に対する算出対象方向での更新後の偏微分値を算出せず、スライス選択傾斜磁場方向以外の傾斜磁場方向の磁化成分に対する算出対象方向での更新後の偏微分値を算出する。例えば、変形例２に係る処理回路１１は、Ａ／Ｄ変換器（ＡＤＣ）の予測出力値を算出する場合、∂ｐＭｚを算出せず、ＩＱ信号に対応する∂ｐＭｘ及び∂ｐＭｙのみを算出すればよい。∂ｐＭｚを算出する計算時間を削減し、Ａ／Ｄ変換器（ＡＤＣ）の予測出力値の算出時間を短縮することが可能になる。
【0129】
なお、処理回路１１は、空間偏微分を用いた連続系のＡＤＣ計算を実行してもよい。これは拡散項を用いた非特許文献２と同様の計算により実行可能である。
【0130】
以上説明した少なくとも１つの実施形態によれば、磁気共鳴数値シミュレーションの計算効率を向上することができる。
【0131】
上記説明において用いた「プロセッサ」という文言は、例えば、ＣＰＵ、ＧＰＵ、或いは、特定用途向け集積回路（Application Specific Integrated Circuit：ＡＳＩＣ）、プログラマブル論理デバイス（例えば、単純プログラマブル論理デバイス（Simple Programmable Logic Device：ＳＰＬＤ）、複合プログラマブル論理デバイス（Complex Programmable Logic Device：ＣＰＬＤ）、及びフィールドプログラマブルゲートアレイ（Field Programmable Gate Array：ＦＰＧＡ））等の回路を意味する。プロセッサは記憶回路に保存されたプログラムを読み出し実行することで機能を実現する。なお、記憶回路にプログラムを保存する代わりに、プロセッサの回路内にプログラムを直接組み込むよう構成しても構わない。この場合、プロセッサは回路内に組み込まれたプログラムを読み出し実行することで機能を実現する。一方、プロセッサが例えばＡＳＩＣである場合、プログラムが記憶回路に保存される代わりに、当該機能がプロセッサの回路内に論理回路として直接組み込まれる。なお、本実施形態の各プロセッサは、プロセッサごとに単一の回路として構成される場合に限らず、複数の独立した回路を組み合わせて１つのプロセッサとして構成し、その機能を実現するようにしてもよい。さらに、図１及び図９における複数の構成要素を１つのプロセッサへ統合してその機能を実現するようにしてもよい。
【0132】
いくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更、実施形態同士の組み合わせを行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれると同様に、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。
【符号の説明】
【0133】
１ 磁気共鳴数値シミュレーション装置
１１ 処理回路
１２ 通信機器
１３ 表示機器
１４ 入力機器
１５ 記憶装置
４１ 非印加期間
４２ 印加期間
１１１ 入力機能
１１２ 数値計算機能
１１３ 出力制御機能
１１４ 判定機能
１１５ 行列算出機能
１１ 処理回路
１２ 通信機器
１３ 表示機器
１４ 入力機器
１５ 記憶装置
４１ 非印加期間
４２ 印加期間
１１１ 入力機能
１１２ 数値計算機能
１１３ 出力制御機能
１１４ 判定機能
１１５ 行列算出機能
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
