JP6548615B2 - 磁場歪み算出装置、方法およびプログラム - Google Patents

Description

本発明は、ＭＲＩ（Magnetic Resonance Imaging）装置により取得された３次元画像に含まれる磁場歪みを算出する磁場歪み算出装置、方法およびプログラムに関するものである。

近年、ＣＴ（Computed Tomography）装置およびＭＲＩ装置等の医療機器の進歩により、質の高い高解像度の３次元画像が画像診断に用いられるようになってきている。ここで、ＣＴ装置は、短時間で被検体の体内を撮影することができ、疾患の早期発見および手術計画の立案に用いられることが多い。一方、ＭＲＩ装置は、ＣＴ装置よりも撮影時間を要するものの、被検体の被曝がないというメリットを有する。ＭＲＩ装置は、被検体である人体に傾斜磁場をかけ、人体の組織を構成する原子核スピンが発生する核磁気共鳴信号を計測し、人体内部の頭部および腹部等の形態を３次元的に画像化する装置である。

しかしながら、ＭＲＩ装置により取得されるＭＲＩ画像には、静的磁場における不均一性および傾斜磁場における不完全性による、装置に起因する幾何学的歪み、並びに被検体に起因する幾何学的歪み（以下、磁場歪みと称する）が発生する。このような磁場歪みは、ある程度は許容されるものの、同一被検体についての経過観察を行うための撮影時期が異なる２つのＭＲＩ画像に磁場歪みが含まれると、経過を正確に判断することができない。とくに、被検体がアルツハイマー病の患者の場合、正常者の脳全体の萎縮率が年１％未満であるのに対して、年に１〜３％である。このため、アルツハイマー病の経過観察においては、前回の診断時に取得したＭＲＩ画像と最新のＭＲＩ画像とを比較して、脳のどの部分がどの程度萎縮しているかを正確に認識する必要がある。しかしながら、ＭＲＩ画像に磁場歪みが含まれていると、求められた脳の萎縮が病気の経過に起因するものであるのか、磁場歪みに起因するものであるのかが判別がつかないこととなる。このため、磁気歪みを補正するための各種手法が提案されている。

例えば、傾斜磁場発生装置を構成するコイルのパターンに基づいて空間内の撮像断面における傾斜磁場強度を算出し、これを用いて磁場歪みを算出してＭＲＩ画像の磁場歪みを補正する手法が提案されている（特許文献１参照）。また、歪みを測定するためのファントムを用いて磁場歪みを算出し、算出した磁場歪みを用いて被検体を撮影することにより取得したＭＲＩ画像を補正する手法も提案されている（特許文献２，３参照）。これらの手法を用いることにより、ＭＲＩ画像に含まれる磁場歪みを補正することができる。

一方、上述したように、同一被検体の経過観察を行うためには、撮影時期が異なる２つの画像の位置合わせを精度よく行って、前回からの変更箇所および変更量を検出する必要がある。このため、一般的な画像の位置合わせの手法として、２つの画像の剛体位置合わせを行って画像間の大域的な変形を表す変形ベクトル場を第１の変形情報として取得し、第１の変形情報を初期変形情報として、２つの画像の非剛体位置合わせを行うことにより、２つの画像の局所的な変形を表す変形ベクトル場を取得する手法が提案されている（特許文献４参照）。

特開２０１２−１６１３５４号公報 特開平１１−９７０８号公報 特表２００８−５２１４７１号公報 特開２０１４−１０８３４９号公報

一方、ＭＲＩ装置においては、機器の経年劣化により磁場歪みが悪化する場合がある。このため、特許文献１〜３に記載された手法においては、例えば半年または１年等の一定期間ごとに磁場歪みを計測するキャリブレーションを必要がある。しかしながら、一定期間ごとに磁場歪みを計測することは非常に手間がかかる。また、キャリブレーション直後においては、ＭＲＩ画像の磁場歪みを精度よく補正することができるが、キャリブレーションを行ってからの経過時間が長いと、実際の磁場歪みと計測された磁場歪みとが異なる可能性が非常に高いため、磁場歪みを精度よく補正することができない。また、特許文献４に記載された手法は、ＭＲＩ画像の磁気歪みの補正を行うことを意図したものではなく、剛体位置合わせおよび非剛体位置合わせを組み合わせてはいるものの、ＭＲＩ画像の磁場歪みを補正した位置合わせを行うことはできない。

本発明は、上記事情に鑑みなされたものであり、キャリブレーションを行うことなく、ＭＲＩ画像に含まれる磁場歪みを補正することを目的とする。

本発明による磁場歪み算出装置は、被検体の対象部位をＭＲＩ装置により撮影にすることにより取得された、撮影時期が異なる第１の３次元画像および第２の３次元画像を取得する画像取得手段と、
第１の３次元画像および第２の３次元画像のそれぞれから、互いに対応する第１の形状不変領域および第２の形状不変領域を抽出する領域抽出手段と、
第１の形状不変領域および第２の形状不変領域を剛体位置合わせして、第１の形状不変領域上の各点の第２の形状不変領域上の対応する各点に対する第１の変形ベクトルを取得する第１の位置合わせ手段と、
第１の３次元画像と第２の３次元画像とを、第１の形状不変領域および第２の形状不変領域において非剛体位置合わせして、第１の３次元画像の第１の形状不変領域上の各点の、第２の３次元画像の第２の形状不変領域上の対応する各点に対する第２の変形ベクトルを取得する第２の位置合わせ手段と、
第１の変形ベクトルおよび第２の変形ベクトルに基づいて、第１の３次元画像上の各点における、第１の３次元画像および第２の３次元画像間の相対的な磁場歪みを表す磁場歪みベクトルを算出する磁場歪みベクトル算出手段とを備えたことを特徴とするものである。

「剛体位置合わせ」とは、位置合わせの対象が変形しないと仮定した場合における位置合わせである。例えば３次元画像全体で剛体位置合わせを行うことにより、第１の３次元画像上の各点における第２の３次元画像に対する平行移動および回転を表すベクトルを第１の変形ベクトルとして取得することができる。すなわち、剛体位置合わせは、３次元画像間における局部の非線形な変形の位置合わせは含まない。

「非剛体位置合わせ」とは、位置合わせの対象が変形する場合における位置合わせである。非剛体位置合わせにより、３次元画像間の非線形な変形（例えば画像の局所的なねじれ）も含めて位置合わせすることが可能である。なお、非剛体位置合わせする領域は、第１の形状不変領域および第２の形状不変領域のみであってもよく、第１の形状不変領域および第２の形状不変領域を含む、第１の形状不変領域および第２の形状不変領域内の領域であってもよい、
「各点」とは、３次元画像上の全ての画素位置（ボクセル位置）であってもよく、あらかじめ定められた間隔で間引いた画素位置であってもよい。

なお、本発明による磁場歪み算出装置においては、磁場歪みベクトル算出手段は、第１の変形ベクトルに基づいて第１の形状不変領域における第２の変形ベクトルを補正することにより取得される第３の変形ベクトルを第１の形状不変領域における磁場歪みベクトルとして算出し、第１の３次元画像における第１の形状不変領域以外の他の領域においては、第３の変形ベクトルに基づいて磁場歪みベクトルを算出するものであってもよい。

また、本発明による磁場歪み算出装置においては、磁場歪みベクトルに基づいて、第１の３次元画像に含まれる磁場歪みを除去し、磁場歪みが除去された第１の３次元画像および第２の３次元画像を非剛体位置合わせして、第１の３次元画像に含まれる対象部位の各点の、第２の３次元画像に含まれる対象部位の対応する各点に対する第４の変形ベクトルを算出する第３の位置合わせ手段をさらに備えるものとしてもよい。

この場合、第４の変形ベクトルに基づいて、対象部位の容積変化量を算出する変化量算出手段をさらに備えるものとしてもよい。

本発明による磁場歪み算出装置においては、対象部位が脳であってもよく、この場合、第１の形状不変領域および第２の形状不変領域が頭蓋骨であってもよい。

本発明による磁場歪み算出方法は、被検体の対象部位をＭＲＩ装置により撮影にすることにより取得された、撮影時期が異なる第１の３次元画像および第２の３次元画像を取得し、
第１の３次元画像および第２の３次元画像のそれぞれから、互いに対応する第１の形状不変領域および第２の形状不変領域を抽出し、
第１の形状不変領域および第２の形状不変領域を剛体位置合わせして、第１の形状不変領域上の各点の第２の形状不変領域上の対応する各点に対する第１の変形ベクトルを取得し、
第１の３次元画像と第２の３次元画像とを、第１の形状不変領域および第２の形状不変領域において非剛体位置合わせして、第１の３次元画像の第１の形状不変領域上の各点の、第２の３次元画像の第２の形状不変領域上の対応する各点に対する第２の変形ベクトルを取得し、
第１の変形ベクトルおよび第２の変形ベクトルに基づいて、第１の３次元画像上の各点における、第１の３次元画像および第２の３次元画像間の相対的な磁場歪みを表す磁場歪みベクトルを算出することを特徴とするものである。

なお、本発明による磁場歪み算出方法をコンピュータに実行させるためのプログラムとして提供してもよい。

本発明によれば、第１の３次元画像および第２の３次元画像のそれぞれから、互いに対応する第１の形状不変領域および第２の形状不変領域が抽出され、第１の形状不変領域および第２の形状不変領域とが剛体位置合わせされて、第１の形状不変領域の第２の形状不変領域に対する第１の変形ベクトルが取得される。また、第１の３次元画像と第２の３次元画像とが、第１の形状不変領域および第２の形状不変領域において非剛体位置合わせされて、第１の３次元画像上の各点における第２の３次元画像に対する第２の変形ベクトルが取得される。次いで、第１の変形ベクトルおよび第２の変形ベクトルに基づいて、第１の３次元画像上の各点における、第１の３次元画像および第２の３次元画像間の相対的な磁場歪みを表す磁場歪みベクトルが算出される。

ここで、被検体内の第１および第２の形状不変領域は、実際には変形しない領域であるが、第１および第２の３次元画像の取得時の磁場歪みの影響により、第１の３次元画像に含まれる第１の形状不変領域は、第２の３次元画像に含まれる第２の形状不変領域に対して変形している。本発明において取得される第１の変形ベクトルは、第１の形状不変領域の第２の形状不変領域に対する平行移動および回転による位置ずれを表すものとなる。一方、第２の変形ベクトルは、第１の形状不変領域の第２の形状不変領域に対する平行移動および回転による位置ずれ、並びに第１および第２の３次元画像間の相対的な磁場歪みによる変形を含んでいる。

このため、第１の変形ベクトルおよび第２の変形ベクトルに基づくことにより、ファントム等を用いて定期的に磁場歪みを測定するキャリブレーションを行わなくても、第１および第２の３次元画像間における相対的な磁場歪みによる変形、すなわち磁気歪みベクトルを算出することができる。

本発明の第１の実施形態による磁場歪み算出装置を適用した、診断支援システムの概要を示すハードウェア構成図 本実施形態による磁場歪み算出装置の構成を示す概略ブロック図 形状不変領域の抽出を説明するための図 第１の変形ベクトルの取得を説明するための図 第１および第２の形状不変領域Ａ１，Ａ２上の対応する画素位置Ｐ１，Ｐ２と第１の変形ベクトルＶ１との関係を示す図 第２の変形ベクトルの取得を説明するための図 第１および第２の形状不変領域Ａ１，Ａ２上の対応する画素位置Ｐ１，Ｐ２と第２の変形ベクトルＶ２との関係を示す図 第１および第２の３次元画像と第２の変形ベクトルとの関係を説明するための図 第３の変形ベクトルの取得を説明するための図 脳の容積変化量の算出を説明するための図 脳の容積変化量の算出を説明するための図 脳の表面を表す脳画像に容積変化量を可視化した状態を示す図 脳の断面を表す脳画像に容積変化量を可視化した状態を示す図 脳の表面を表す脳画像において、脳区域毎に異常領域を可視化した状態を示す図 本実施形態において行われる処理を示すフローチャート

以下、図面を参照して本発明の実施形態について説明する。図１は、本発明の第１の実施形態による磁場歪み算出装置を適用した、診断支援システムの概要を示すハードウェア構成図である。図１に示すように、診断支援システムでは、本実施形態による磁場歪み算出装置１、３次元画像撮影装置２、および画像保管サーバ３が、ネットワーク４を経由して通信可能な状態で接続されている。そして、診断支援システムにおいては、被検体の診断の対象部位についての比較診断のために、撮影時期が異なる２つの３次元画像の比較を磁場歪み算出装置１において行うものである。

３次元画像撮影装置２は、被検体の診断対象となる部位を撮影することにより、その部位を表す３次元画像を生成する装置であり、本実施形態においては、ＭＲＩ装置である。この３次元画像撮影装置２により生成された３次元画像は画像保管サーバ３に送信され、保存される。なお、本実施形態においては、被検体の診断対象部位は脳であり、３次元画像撮影装置２において、被検体の頭部についての３次元画像が生成されるものとする。

ここで、ＭＲＩ装置により取得される３次元画像には磁場歪みが含まれる。磁場歪みは、ある程度は許容されるものの、撮影時期が異なる２つの３次元画像に磁場歪みが含まれると、対象部位の病状の経過を正確に判断することができない。とくに、被検体がアルツハイマー病の患者の場合、正常者の脳全体の萎縮率が年１％未満であるのに対して、年に１〜３％である。このため、アルツハイマー病の経過観察においては、前回の診断時に取得した３次元画像と最新の３次元画像とを比較して、脳のどの部分がどの程度萎縮しているかを正確に認識する必要がある。しかしながら、３次元画像に磁場歪みが含まれていると、求められた脳の萎縮が病気の経過に起因するものであるのか、磁場歪みに起因するものであるのかが判別がつかないこととなる。

本実施形態においては、最新の３次元画像および過去の３次元画像間の相対的な磁場歪みを表す磁場歪みベクトルを算出し、磁場歪みベクトルを用いて、磁場歪みの影響を受けることなく、脳の容積変化量、すなわち脳の萎縮の程度を診断できるようにするものである。

画像保管サーバ３は、各種データを保存して管理するコンピュータであり、大容量外部記憶装置およびデータベース管理用ソフトウェアを備えている。画像保管サーバ３は、有線あるいは無線のネットワーク４を介して他の装置と通信を行い、画像データ等を送受信する。具体的には３次元画像撮影装置２で生成された３次元画像等の画像データをネットワーク経由で取得し、大容量外部記憶装置等の記録媒体に保存して管理する。なお、画像データの格納形式およびネットワーク４経由での各装置間の通信は、ＤＩＣＯＭ（Digital Imaging and COmmunication in Medicine）等のプロトコルに基づいている。また、本実施形態においては、同一被検体についての撮影時期が異なる頭部の３次元画像が、画像保管サーバ３に保管されているものとする。

磁場歪み算出装置１は、１台のコンピュータに、本発明の磁場歪み算出プログラムをインストールしたものである。コンピュータは、診断を行う医師が直接操作するワークステーションまたはパーソナルコンピュータでもよいし、それらとネットワークを介して接続されたサーバコンピュータでもよい。磁場歪み算出プログラムは、ＤＶＤ（Digital Versatile Disc）あるいはＣＤ−ＲＯＭ（Compact Disk Read Only Memory）等の記録媒体に記録されて配布され、その記録媒体からコンピュータにインストールされる。または、ネットワークに接続されたサーバコンピュータの記憶装置、もしくはネットワークストレージに、外部からアクセス可能な状態で記憶され、要求に応じて医師が使用するコンピュータにダウンロードされ、インストールされる。

図２は、コンピュータに磁場歪み算出プログラムをインストールすることにより実現される磁場歪み算出装置の概略構成を示す図である。図２に示すように、磁場歪み算出装置１は、標準的なワークステーションの構成として、ＣＰＵ（Central Processing Unit）１１、メモリ１２およびストレージ１３を備えている。また、磁場歪み算出装置１には、ディスプレイ１４および、マウス等の入力部１５が接続されている。

ストレージ１３には、ネットワーク４を経由して画像保管サーバ３から取得した、同一被検体についての最新および前回の診断時における過去の３次元画像、並びに処理に必要な情報を含む各種情報が記憶されている。なお、本実施形態においては、同一被検体についての頭部を対象部位とする最新の３次元画像Ｇ１（第１の３次元画像）、および例えば前回の診断時に取得した過去の３次元画像Ｇ２（第２の３次元画像）が記憶されているものとする。

また、メモリ１２には、磁場歪み算出プログラムが記憶されている。磁場歪み算出プログラムは、ＣＰＵ１１に実行させる処理として、３次元画像撮影装置２が取得した、第１および第２の３次元画像Ｇ１，Ｇ２を取得する画像取得処理、第１の３次元画像Ｇ１および第２の３次元画像Ｇ２のそれぞれから、互いに対応する第１の形状不変領域および第２の形状不変領域を抽出する領域抽出処理、第１の形状不変領域および第２の形状不変領域を剛体位置合わせして、第１の形状不変領域上の各点の第２の形状不変領域上の対応する各点に対する第１の変形ベクトルを取得する第１の位置合わせ処理、第１の３次元画像Ｇ１と第２の３次元画像Ｇ２とを、第１の形状不変領域および第２の形状不変領域において非剛体位置合わせして、第１の３次元画像Ｇ１の第１の形状不変領域上の各点の、第２の３次元画像Ｇ２の第２の形状不変領域上の対応する各点に対する第２の変形ベクトルを取得する第２の位置合わせ処理、第１の変形ベクトルおよび第２の変形ベクトルに基づいて、第１の３次元画像Ｇ１上の各点における、第１および第２の３次元画像Ｇ１，Ｇ２間の相対的な磁場歪みを表す磁場歪みベクトルを算出する磁場歪みベクトル算出処理、磁場歪みベクトルに基づいて、第１の３次元画像Ｇ１に含まれる磁場歪みを除去し、磁場歪みが除去された第１の３次元画像Ｇ１１および第２の３次元画像Ｇ２を位置合わせして、第１の３次元画像Ｇ１に含まれる対象部位である脳の各点の、第２の３次元画像Ｇ２に含まれる脳の対応する各点に対する第４の変形ベクトルを算出する第３の位置合わせ処理、並びに、第４の変形ベクトルに基づいて、脳の容積変化量を算出する変化量算出処理を規定する。

そして、ＣＰＵ１１がプログラムに従いこれらの処理を実行することで、コンピュータは、画像取得部２１、領域抽出部２２、第１の位置合わせ部２３、第２の位置合わせ部２４、磁場歪みベクトル算出部２５、第３の位置合わせ部２６、および変化量算出部２７として機能する。なお、磁場歪み算出装置１は、画像取得処理、領域抽出処理、第１の位置合わせ処理、第２の位置合わせ処理、磁場歪みベクトル算出処理、第３の位置合わせ処理、および変化量算出処理をそれぞれ行う複数のプロセッサまたは処理回路を備えるものであってもよい。

画像取得部２１は、対象部位である脳を含む頭部の第１および第２の３次元画像Ｇ１，Ｇ２を画像保管サーバ３から取得する。なお、第１および第２の３次元画像Ｇ１，Ｇ２が既にストレージ１３に記憶されている場合には、画像取得部２１は、ストレージ１３から第１および第２の３次元画像Ｇ１，Ｇ２を取得するようにしてもよい。

領域抽出部２２は、第１の３次元画像Ｇ１および第２の３次元画像Ｇ２のそれぞれから、互いに対応する第１の形状不変領域および第２の形状不変領域を抽出する。本実施形態においては、領域抽出部２２は、第１の３次元画像Ｇ１および第２の３次元画像Ｇ２のそれぞれから、頭蓋骨の領域を第１の形状不変領域Ａ１および第２の形状不変領域Ａ２として抽出する。ここで、頭蓋骨は、第１および第２の３次元画像Ｇ１，Ｇ２において、比較的高輝度の領域として含まれる。領域抽出部２２は、第１および第２の３次元画像Ｇ１，Ｇ２における骨と見なせる比較的高輝度の画素値（ボクセル値）を有する領域を、しきい値処理により第１および第２の形状不変領域Ａ１，Ａ２として抽出する。図３は第１および第２の３次元画像Ｇ１，Ｇ２からの第１および第２の形状不変領域Ａ１，Ａ２の抽出を説明するための図である。図３に示すように、領域抽出部２２は、第１の３次元画像Ｇ１に含まれる頭部の組織から頭蓋骨のみを第１の形状不変領域Ａ１として抽出する。同様に、領域抽出部２２は、第２の３次元画像Ｇ２についても、同様に第２の形状不変領域Ａ２を抽出する。なお、第１および第２の形状不変領域Ａ１，Ａ２の抽出方法はこれに限定されるものではなく、ＭＲＩ装置により取得した３次元画像から骨領域を抽出するように学習された判別器を用いる手法等、任意の手法を用いることができる。また、頭蓋骨には、動きがある下顎骨は含まないものとする。

第１の位置合わせ部２３は、第１の形状不変領域Ａ１および第２の形状不変領域Ａ２を剛体位置合わせして、第１の形状不変領域Ａ１上の各点の第２の形状不変領域Ａ２上の対応する各点に対する第１の変形ベクトルＶ１を取得する第１の位置合わせ処理を行う。具体的には、第１の位置合わせ部２３は、最小二乗法を用いて、第１の形状不変領域Ａ１上の各点と、これらの各点と対応する第２の形状不変領域Ａ２上の各点との距離の二乗和が最小となるように、第１の形状不変領域Ａ１および第２の形状不変領域Ａ２の対応する各点間の平行移動成分および回転成分を算出する。なお、各点は、第１および第２の形状不変領域Ａ１，Ａ２における全ての画素位置であってもよく、適宜間引いた画素位置であってもよい。また、平行移動成分および回転成分の算出は、最小二乗法に限定されるものではなく、任意の手法を用いることができる。

図４は第１の変形ベクトルＶ１の取得を説明するための図である。図４に示すように、第１の位置合わせ部２３は、第１の形状不変領域Ａ１上の各点において、平行移動成分および回転成分を変形ベクトルＶｒ１として算出する。なお、図４においては、説明のために、第１の形状不変領域Ａ１上の適宜間引いた位置において取得された変形ベクトルＶｒ１を示している。なお、図４および以降の説明においては、変形ベクトルの取得を２次元画像において示しているが、実際には３次元画像において３次元の変形ベクトルが取得される。

そして、第１の位置合わせ部２３は、図４に示すように、各画素の平行移動成分および回転成分を表す変形ベクトルＶｒ１の平均値を算出し、第１の変形ベクトルＶ１を取得する。これにより、第１の変形ベクトルＶ１は、第１の形状不変領域Ａ１の第２の形状不変領域Ａ２に対する平行移動および回転による位置ずれを表すものとなる。

図５は、第１および第２の形状不変領域Ａ１，Ａ２上の対応する画素位置Ｐ１，Ｐ２と第１の変形ベクトルＶ１との関係を示す図である。図５に示すように、第１の形状不変領域Ａ１における画素位置Ｐ１から、第２の形状不変領域Ａ２における画素位置Ｐ２への座標変換は、第１の変形ベクトルＶ１を用いて下記の式（１）により表される。なお、図５の下の図においては、第１の形状不変領域Ａ１を破線で、第２の形状不変領域Ａ２を実線で示している。

Ｐ２＝Ｐ１＋Ｖ１ （１）
第２の位置合わせ部２４は、第１の３次元画像Ｇ１と第２の３次元画像Ｇ２とを、第１の形状不変領域Ａ１および第２の形状不変領域Ａ２において非剛体位置合わせして、第１の３次元画像Ｇ１の第１の形状不変領域Ａ１上の各点の、第２の３次元画像Ｇ２の第２の形状不変領域Ａ２上の対応する各点に対する第２の変形ベクトルＶ２を取得する、第２の位置合わせを行う。本実施形態においては、第１の形状不変領域Ａ１と第２の形状不変領域Ａ２とを非剛体位置合わせして、第２の変形ベクトルＶ２を取得する。なお、第１の３次元画像Ｇ１の第１の形状不変領域Ａ１を含む第１の形状不変領域Ａ１内の領域と、第２の３次元画像Ｇ２における第２の形状不変領域Ａ２を含む第２の形状不変領域Ａ２内の領域とを非剛体位置合わせし、これにより取得される変形ベクトルのうち、第１の形状不変領域Ａ１と第２の形状不変領域Ａ２とに存在するもののみを第２の変形ベクトルとして取得してもよい。

なお、本実施形態における非剛体位置合わせとは、第１の形状不変領域Ａ１上の各点を各点の変形量に基づいて移動させて、第１および第２の形状不変領域Ａ１，Ａ２の類似度を判定する所定の関数を最大化または最小化することにより、第１および第２の形状不変領域Ａ１，Ａ２を互いに一致させるための、第１の形状不変領域Ａ１の各点の変形量を、第２の変形ベクトルＶ２として算出して取得する方法である。非剛体位置合わせとしては、例えば、Rueckert D Sonoda LI，Hayesc，Ｅt al.、「Nonrigid Registration Using Free-Form Deformations：application to breast MR Images」、IEEE transactions on Medical Imaging、1999年、vol.18，No.8，pp.712-721等種々の周知の手法を適用することができる。

図６は第２の変形ベクトルＶ２の取得を説明するための図である。図６に示すように、第２の位置合わせ部２４は、第１の形状不変領域Ａ１上の各点において、第２の変形ベクトルＶ２を取得する。なお、図６においては、説明のために、第１の形状不変領域Ａ１上の適宜間引いた位置において取得された第２の変形ベクトルＶ２を示している。

図７は、第１および第２の形状不変領域Ａ１，Ａ２上の対応する画素位置Ｐ１，Ｐ２と第２の変形ベクトルＶ２との関係を示す図である。図７に示すように、第１の形状不変領域Ａ１における画素位置Ｐ１から、第２の形状不変領域Ａ２における画素位置Ｐ２への座標変換は、第２の変形ベクトルＶ２を用いて下記の式（２）により表される。なお、図７の下の図においても、第１の形状不変領域Ａ１を破線で、第２の形状不変領域Ａ２を実線で示している。

Ｐ２＝Ｐ１＋Ｖ２ （２）
図８は、第１および第２の３次元画像Ｇ１，Ｇ２と第２の変形ベクトルＶ２との関係を説明するための図である。第１の３次元画像Ｇ１の第１の形状不変領域Ａ１を第２の変形ベクトルＶ２により変形すると、第２の３次元画像Ｇ２の第２の形状不変領域Ａ２が得られる。

一方、本実施形態においては、第１および第２の３次元画像Ｇ１，Ｇ２は、ＭＲＩ装置により取得されたものであるため、磁場歪みの影響により、各３次元画像Ｇ１，Ｇ２に含まれる被検体の頭蓋骨の領域は、実際の頭蓋骨からみて歪んだものとなっている。また、第１の３次元画像Ｇ１と第２の３次元画像Ｇ２とでは、撮影時期の相違により、磁場歪みの程度が異なるため、頭蓋骨の変形の程度が異なる。また、アルツハイマー病の患者の病気が進行して脳が萎縮するため、第１の３次元画像Ｇ１と第２の３次元画像Ｇ２とでは、撮影時期の相違による脳の萎縮により、脳が変形している。脳の萎縮に起因する変形を求め、アルツハイマー病の診断を精度よく行うためには、第１の３次元画像Ｇ１上の各点における、第１および第２の３次元画像Ｇ１，Ｇ２間の相対的な磁場歪みを除去した上で、第１および第２の３次元画像Ｇ１，Ｇ２の位置合わせを行う必要がある。

ここで、第１および第２の形状不変領域Ａ１，Ａ２は、下顎骨を除く頭蓋骨である。アルツハイマー病の診断を行うのは年配の患者が多く、時間の経過によって頭蓋骨が成長して変形することはほぼない。このため、第１および第２の形状不変領域Ａ１，Ａ２から算出した第２の変形ベクトルＶ２は、頭蓋骨自体の変形は含まないが、第１の形状不変領域Ａ１の第２の形状不変領域Ａ２に対する平行移動および回転による位置ずれＰ０、並びに第１および第２の３次元画像Ｇ１，Ｇ２間の相対的な磁場歪みによる変形Ｐｍを含んでいる。

磁場歪みベクトル算出部２５は、第１の変形ベクトルＶ１および第２の変形ベクトルＶ２に基づいて、第１の３次元画像Ｇ１上の各点における、第１および第２の３次元画像Ｇ１，Ｇ２間の相対的な磁場歪みを表す磁場歪みベクトルＶｍを算出する。具体的には、第１の変形ベクトルＶ１に基づいて第１の形状不変領域Ａ１における第２の変形ベクトルＶ２を補正することにより、第３の変形ベクトルＶ３を第１の形状不変領域Ａ１における磁場歪みベクトルＶｍとして算出する。ここで、第１の変形ベクトルＶ１は、第１の形状不変領域Ａ１の第２の形状不変領域Ａ２に対する平行移動および回転による位置ずれを表す。このため、第３の変形ベクトルＶ３の算出は、第２の変形ベクトルＶ２から第１の変形ベクトルＶ１を減算することにより行えばよい。

図９は第３の変形ベクトルＶ３の取得を説明するための図である。図９に示すように、第１の位置合わせ部２３は、第１の形状不変領域Ａ１上の各点において、第３の変形ベクトルＶ３を磁場歪みベクトルＶｍとして算出する。さらに、磁場歪みベクトル算出部２５は、第３の変形ベクトルＶ３に基づいて、第１の３次元画像Ｇ１における第１の形状不変領域Ａ１以外の他の領域Ａ３における磁場歪みベクトルＶｍを算出する。具体的には、アフィン変換モデル、シンプレートスプラインモデル等の公知のモデルへの関数フィッティングを行うことにより、第１の形状不変領域Ａ１上の第３の変形ベクトルＶ３を空間的に補間して、他の領域Ａ３における磁場歪みベクトルＶｍを算出する。なお、他の領域Ａ３における磁場歪みベクトルＶｍの算出の手法は、これらに限定されるものではなく、任意の補間の手法を用いることができる。これにより、図９に示すように、第１の形状不変領域Ａ１および他の領域Ａ３における磁場歪みベクトルＶｍが算出される。なお、このようにして算出した磁場歪みベクトルＶｍは、図８における第１および第２の３次元画像Ｇ１，Ｇ２間の相対的な磁場歪みによる変形Ｐｍと一致する。

第３の位置合わせ部２６は、磁場歪みベクトルＶｍに基づいて、第１の３次元画像Ｇ１に含まれる磁場歪みを除去する。具体的には、図８に示すように、磁場歪みベクトルＶｍに基づいて第１の３次元画像Ｇ１を変形することにより、磁場歪みが除去された第１の３次元画像Ｇ１１を生成する。そして、第３の位置合わせ部２６は、磁場歪みが除去された第１の３次元画像Ｇ１１を第２の３次元画像Ｇ２と非剛体位置合わせして、第１の３次元画像Ｇ１１に含まれる対象部位である脳の各点の、第２の３次元画像Ｇ２に含まれる脳の対応する各点に対する第４の変形ベクトルＶ４を算出する第３の位置合わせ処理を行う。非剛体位置合わせは、上記第２の位置合わせ部２４と同様に行えばよい。なお、磁場歪みの除去および第４の変形ベクトルＶ４の算出は、第１の３次元画像Ｇ１１における頭蓋骨および頭蓋骨内部の領域についてのみ行えばよいが、第１の３次元画像Ｇ１１の全体において行ってもよい。ここで、第４の変形ベクトルＶ４は、第１および第２の３次元画像Ｇ１，Ｇ２の撮影時期の相違による、脳の萎縮による変形のみを表すものとなる。したがって、図８に示すように、磁場歪みが除去された第１の３次元画像Ｇ１１を第４の変形ベクトルＶ４により変形すると、第２の３次元画像Ｇ２が得られることとなる。

なお、磁場歪みが除去された第１の３次元画像Ｇ１１内の第１の形状不変領域Ａ１と、第２の３次元画像Ｇ２内の第２の形状不変領域Ａ２とは、剛体位置合わせにより表される平行移動成分および回転成分により、一致させることができる。

変化量算出部２７は、第４の変形ベクトルＶ４に基づいて、脳の容積変化量を算出する。以下、脳の容積変化量の算出について説明する。変化量算出部２７は、まず、第１の３次元画像Ｇ１から脳領域を抽出する。そして、脳領域内の各画素位置において、容積変化量を算出する。図１０および図１１は脳の容積変化量の算出を説明するための図である。なお、脳領域内のある画素位置（ｘ，ｙ，ｚ）における第４の変形ベクトルＶ４をＶ４（ｘ，ｙ，ｚ）とする。まず、画素位置（ｘ，ｙ，ｚ）を、ｘ軸、ｙ軸およびｚ軸の方向にそれぞれｄｘ、ｄｙ、ｄｚ変位させ、これにより得られる（ｘ，ｙ，ｚ）、（ｘ＋ｄｘ，ｙ，ｚ）、（ｘ，ｙ＋ｄｙ，ｚ）、（ｘ，ｙ，ｚ＋ｄｚ）、（ｘ＋ｄｘ，ｙ＋ｄｙ，ｚ）、（ｘ＋ｄｘ，ｙ，ｚ＋ｄｚ）、（ｘ，ｙ＋ｄｙ，ｚ＋ｄｚ）、（ｘ＋ｄｘ，ｙ＋ｄｙ，ｚ＋ｄｚ）の８点を頂点とする６面体を考える。この６面体は直方体であり、その体積ＶＯＬ１はｄｘ×ｄｙ×ｄｚにより算出される。

一方、上記の８点を第４の変形ベクトルｖ４により変形すると、それぞれ｛（ｘ，ｙ，ｚ）＋Ｖ４（ｘ，ｙ，ｚ）｝、｛（ｘ＋ｄｘ，ｙ，ｚ）＋Ｖ４（ｘ＋ｄｘ，ｙ，ｚ）｝、｛（ｘ，ｙ＋ｄｙ，ｚ）＋Ｖ４（ｘ，ｙ＋ｄｙ，ｚ）｝、｛（ｘ，ｙ，ｚ＋ｄｚ）＋Ｖ４（ｘ，ｙ，ｚ＋ｄｚ）｝、｛（ｘ＋ｄｘ，ｙ＋ｄｙ，ｚ）＋Ｖ４（ｘ＋ｄｘ，ｙ＋ｄｙ，ｚ）｝、｛（ｘ＋ｄｘ，ｙ，ｚ＋ｄｚ）＋Ｖ４（ｘ＋ｄｘ，ｙ，ｚ＋ｄｚ）｝、｛（ｘ，ｙ＋ｄｙ，ｚ＋ｄｚ）＋Ｖ４（ｘ，ｙ＋ｄｙ，ｚ＋ｄｚ）｝、｛（ｘ＋ｄｚ，ｙ＋ｄｙ，ｚ＋ｄｚ）＋Ｖ４（ｘ＋ｄｚ，ｙ＋ｄｙ，ｚ＋ｄｚ｝となり、図１０に示す直方体は、例えば図１１に示すように変形される。変化量算出部２７は、変形された６面体の体積ＶＯＬ２を算出する。具体的には、変化量算出部２７は、８点の画素位置の座標値の平均値を算出する。平均値となる画素位置は変形された６面体の内部の点となる。そして、平均値となる画素位置を頂点とする６つの四角錐の体積の合計値を、変形された６面体の体積ＶＯＬ２として算出する。

そして、変化量算出部２７は、各画素位置における容積変化量を算出する。なお、容積変化量はＶＯＬ２／ＶＯＬ１−１により算出する。ここで、容積に変化がなければＶＯＬ２／ＶＯＬ１＝１となる。一方、容積が縮小していればＶＯＬ２／ＶＯＬ１は１より小さい値となるため、容積変化量は負の値となる。容積が膨張していればＶＯＬ２／ＶＯＬ１は１より大きい値となるため、容積変化量は正の値となる。変化量算出部２７は、脳内の各画素位置の容積変化量を脳の萎縮率として、第１の３次元画像Ｇ１から生成した脳画像に可視化して、ディスプレイ１４に表示する。以下、容積変化量の可視化について説明する。

図１２は、第１の３次元画像Ｇ１から生成した脳の表面を表す画像に容積変化量を可視化した状態を示す図である。なお、図１２は脳の左側面の脳画像である。図１２に示すように、脳画像において、容積変化量が±３％を超える画素位置、すなわち容積変化量の絶対値が３％を超える画素位置を含む領域（以下、異常領域とする）Ａ１０が、脳の表面とは異なる色、例えば赤色により示されている。なお、図１２においては色を斜線により示している。また、図１２においては、脳の左側面の画像を示しているが、左側面、右側面、後面、前面、正面および下面のいずれか、またはこれらのうちの複数の画像を表示して、異常領域Ａ１０を可視化してもよい。

図１３は、第１の３次元画像Ｇ１から生成した脳の断面を表す脳画像に容積変化量を可視化した状態を示す図である。図１３に示すように、脳のアキシャル断面の脳画像において、容積変化量の絶対値が３％を超える画素位置を含む異常領域Ａ１０が、その周囲とは異なる色、例えば赤色により示されている。なお、図１３においては色を斜線により示している。なお、表示された断面の位置は、操作者の操作により変更が可能である。また、図１３においては、脳のアキシャル断面の画像を示しているが、アキシャル断面、コロナル断面およびサジタル断面のいずれか、またはこれらのうちの複数の断面の画像を表示して、異常領域Ａ１０を可視化してもよい。

なお、図１２および図１３においては、容積変化量の絶対値が３％を超える異常領域Ａ１０を１つの色により可視化しているが、容積変化量に応じて段階的に異なる色により異常領域を可視化してもよい。この場合、例えば、容積変化量の絶対値が３％を超える領域を赤色、２％以上３％未満の領域をオレンジ色、１％以上２％未満の領域を黄色、１％未満の領域を緑色として、容積変化量を可視化してもよい。また、容積変化量に応じて段階的に異なる色により各画素位置を可視化してもよい。また、脳の断面を表す脳画像において、各画素位置を容積変化量に応じて段階的に異なる色により可視化し、さらに断面位置を順次変更するように動画像を生成してもよい。これにより、脳の各位置における容積変化量を動画像上において確認することができる。

一方、脳は、大脳新皮質を解剖学的に複数の脳区域に分割することができる。このため、抽出した脳を複数の脳区域に分割し、脳区域毎に容積変化量を算出してもよい。この場合、各脳区域において算出した各画素位置の容積変化量の平均値、最大値および標準偏差等の代表値を算出し、代表値の絶対値が３％を超える脳区域を異常領域Ａ１０として、第１の３次元画像Ｇ１に可視化してもよい。図１４は脳の表面を表す脳画像において、脳区域毎に容積変化量を可視化した状態を示す図である。なお、図１４においては、脳の表面を表す画像を脳区域に分割しているが、脳の断面を表す脳画像を脳区域に分割して脳区域毎に異常領域を可視化してもよい。

次いで、本実施形態において行われる処理について説明する。図１５は本実施形態において行われる処理を示すフローチャートである。まず、画像取得部２１が、第１および第２の３次元画像Ｇ１，Ｇ２を取得し（ステップＳＴ１）、領域抽出部２２が、第１の３次元画像Ｇ１および第２の３次元画像Ｇ１のそれぞれから、互いに対応する第１の形状不変領域Ａ１および第２の形状不変領域Ａ２を抽出する（ステップＳＴ２）。次いで、第１の位置合わせ部２３が、第１の形状不変領域Ａ１および第２の形状不変領域Ａ２を剛体位置合わせして、第１の形状不変領域Ａ１上の各点の第２の形状不変領域Ａ２上の対応する各点に対する第１の変形ベクトルＶ１を取得する第１の位置合わせ処理を行う（ステップＳＴ３）。

続いて、第２の位置合わせ部２４が、第１の３次元画像Ｇ１と第２の３次元画像Ｇ２とを、第１の形状不変領域Ａ１および第２の形状不変領域Ａ２において非剛体位置合わせして、第２の変形ベクトルＶ２を取得する、第２の位置合わせ処理を行う（ステップＳＴ４）。そして、磁場歪みベクトル算出部２５が、第１の変形ベクトルＶ１および第２の変形ベクトルＶ２に基づいて、第１の３次元画像Ｇ１上の各点における、第１および第２の３次元画像Ｇ１，Ｇ２間の相対的な磁場歪みを表す磁場歪みベクトルＶｍを算出する（ステップＳＴ５）。次いで、第３の位置合わせ部２６が、磁場歪みベクトルＶｍに基づいて、第１の３次元画像Ｇ１に含まれる磁場歪みを除去して、磁場歪みが除去された第１の３次元画像Ｇ１１を生成する（ステップＳＴ６）。そして、第３の位置合わせ部２６が、磁場歪みが除去された第１の３次元画像Ｇ１１を第２の３次元画像Ｇ２と非剛体位置合わせして、第１の３次元画像Ｇ１１に含まれる対象部位である脳の各点の、第２の３次元画像Ｇ２に含まれる脳の対応する各点に対する第４の変形ベクトルＶ４を算出する第３の位置合わせ処理を行う（ステップＳＴ７）。さらに、変化量算出部２７が、第４の変形ベクトルＶ４に基づいて、脳の容積変化量を算出し（ステップＳＴ８）、算出した脳の容積変化量を第１の３次元画像Ｇ１から得られる脳画像上に可視化してディスプレイ１４に表示し（ステップＳＴ９）、処理を終了する。

このように、本実施形態においては、第１の変形ベクトルＶ１および第２の変形ベクトルＶ２に基づいて、第１の３次元画像Ｇ１上の各点における、第１および第２の３次元画像Ｇ１，Ｇ２間の相対的な磁場歪みを表す磁場歪みベクトルＶｍを算出するようにしたものである。このため、ファントム等を用いて定期的に磁場歪みを測定するキャリブレーションを行わなくても、第１および第２の３次元画像Ｇ１，Ｇ２間における相対的な磁場歪みベクトルＶｍを算出することができる。

また、磁場歪みベクトルＶｍに基づいて、第１の３次元画像Ｇ１に含まれる磁場歪みを除去し、磁場歪みが除去された第１の３次元画像Ｇ１１および第２の３次元画像Ｇ２を非剛体位置合わせして、第４の変形ベクトルＶ４を算出することにより、磁場歪みによる影響を除去して、第１および第２の３次元画像Ｇ１，Ｇ２間における対象部位である脳の各位置における変形のみを第４の変形ベクトルＶ４として算出することができる。したがって、第１および第２の３次元画像Ｇ１，Ｇ２間の脳の各位置における経時的な変形を、磁気歪みに影響されることなく、精度よく算出することができる。

また、第４の変形ベクトルＶ４に基づいて、脳の容積変化量を算出することにより、アルツハイマー病の診断に必要な脳の各位置における萎縮率を精度よく算出することができる。

なお、上記実施形態においては、変化量算出部２７において、第４の変形ベクトルＶ４に基づいて脳の容積変化量を算出している。しかしながら、第４の変形ベクトルＶ４を算出することなく、磁場歪みを除去した後の第１の３次元画像Ｇ１１における脳の容積、および第２の３次元画像Ｇ２における脳の容積を各画像のボクセル数から算出し、その容積の差を容積変化量として算出してもよい。

また、上記実施形態においては、最新の３次元画像を過去の３次元画像に位置合わせしているが、過去の３次元画像を第１の３次元画像Ｇ１とし、最新の３次元画像を第２の３次元画像Ｇ２として、過去の３次元画像を最新の３次元画像に位置合わせしてもよい。

以下、本実施形態の作用効果について説明する。

磁場歪みベクトルに基づいて、第１の３次元画像に含まれる磁場歪みを除去し、磁場歪みが除去された第１の３次元画像および第２の３次元画像を非剛体位置合わせして、第１の３次元画像に含まれる対象部位の各点の、第２の３次元画像に含まれる対象部位の対応する各点に対する第４の変形ベクトルを算出することにより、磁場歪みによる影響を除去して、第１および第２の３次元画像間における対象部位の変形のみを第４の変形ベクトルとして算出することができる。したがって、第１および第２の３次元画像間における対象部位の経時的な変形を、磁気歪みに影響されることなく、精度よく算出することができる。

第４の変形ベクトルに基づいて、対象部位の容積変化量を算出することにより、対象部位の容積変化量を精度よく算出することができる。

対象部位を脳とすることにより、第１および第２の３次元画像間の脳の各位置における経時的な変形を算出することができる。したがって、アルツハイマー病の診断に必要な脳の萎縮率を精度よく算出することができる。

頭蓋骨は５０歳を過ぎると成長が止まり、経時による変形はほぼ生じない。このため、アルツハイマー病を発症する高齢者の脳の萎縮率を算出するに際し、形状不変領域を頭蓋骨とすることにより、経時による形状不変領域の変形に影響されることなく、脳の萎縮率を精度よく算出することができる。

１ 磁場歪み算出装置
２ ３次元画像撮影装置
３ 画像保管サーバ
４ ネットワーク
１１ ＣＰＵ
１２ メモリ
１３ ストレージ
１４ ディスプレイ
１５ 入力部
２１ 画像取得部
２２ 領域抽出部
２３ 第１の位置合わせ部
２４ 第２の位置合わせ部
２５ 磁場歪みベクトル算出部
２６ 第３の位置合わせ部
２７ 変化量算出部
Ｇ１，Ｇ２ ３次元画像
Ａ１，Ａ２ 形状不変領域
Ａ１０ 異常領域
Ｖ１，Ｖ２ 変形ベクトル
Ｖｍ 磁場歪みベクトル

Claims (8)

1. 被検体の対象部位をＭＲＩ装置により撮影にすることにより取得された、撮影時期が異なる第１の３次元画像および第２の３次元画像を取得する画像取得手段と、
   前記第１の３次元画像および前記第２の３次元画像のそれぞれから、互いに対応する第１の形状不変領域および第２の形状不変領域を抽出する領域抽出手段と、
   前記第１の形状不変領域および前記第２の形状不変領域を剛体位置合わせして、前記第１の形状不変領域上の各点の前記第２の形状不変領域上の対応する各点に対する第１の変形ベクトルを取得する第１の位置合わせ手段と、
   前記第１の３次元画像と前記第２の３次元画像とを、前記第１の形状不変領域および前記第２の形状不変領域において非剛体位置合わせして、前記第１の３次元画像の前記第１の形状不変領域上の各点の、前記第２の３次元画像の前記第２の形状不変領域上の対応する各点に対する第２の変形ベクトルを取得する第２の位置合わせ手段と、
   前記第１の変形ベクトルおよび前記第２の変形ベクトルに基づいて、前記第１の３次元画像上の各点における、前記第１の３次元画像および前記第２の３次元画像間の相対的な磁場歪みを表す磁場歪みベクトルを算出する磁場歪みベクトル算出手段とを備えたことを特徴とする磁場歪み算出装置。
2. 前記磁場歪みベクトル算出手段は、前記第１の変形ベクトルに基づいて前記第１の形状不変領域における前記第２の変形ベクトルを補正することにより取得される第３の変形ベクトルを前記第１の形状不変領域における前記磁場歪みベクトルとして算出し、前記第１の３次元画像における前記第１の形状不変領域以外の他の領域においては、前記第３の変形ベクトルに基づいて前記磁場歪みベクトルを算出する請求項１記載の磁場歪み算出装置。
3. 前記磁場歪みベクトルに基づいて、前記第１の３次元画像に含まれる磁場歪みを除去し、該磁場歪みが除去された第１の３次元画像および前記第２の３次元画像を非剛体位置合わせして、前記第１の３次元画像に含まれる前記対象部位の各点の、前記第２の３次元画像に含まれる前記対象部位の対応する各点に対する第４の変形ベクトルを算出する第３の位置合わせ手段をさらに備えた請求項１または２記載の磁場歪み算出装置。
4. 前記第４の変形ベクトルに基づいて、前記対象部位の容積変化量を算出する変化量算出手段をさらに備えた請求項３記載の磁場歪み算出装置
5. 前記対象部位が脳である請求項１から４のいずれか１項記載の磁場歪み算出装置。
6. 前記第１の形状不変領域および前記第２の形状不変領域が頭蓋骨である請求項５記載の磁場歪み算出装置。
7. 被検体の対象部位をＭＲＩ装置により撮影にすることにより取得された、撮影時期が異なる第１の３次元画像および第２の３次元画像を取得し、
   前記第１の３次元画像および前記第２の３次元画像のそれぞれから、互いに対応する第１の形状不変領域および第２の形状不変領域を抽出し、
   前記第１の形状不変領域および前記第２の形状不変領域を剛体位置合わせして、前記第１の形状不変領域上の各点の前記第２の形状不変領域上の対応する各点に対する第１の変形ベクトルを取得し、
   前記第１の３次元画像と前記第２の３次元画像とを、前記第１の形状不変領域および前記第２の形状不変領域において非剛体位置合わせして、前記第１の３次元画像の前記第１の形状不変領域上の各点の、前記第２の３次元画像の前記第２の形状不変領域上の対応する各点に対する第２の変形ベクトルを取得し、
   前記第１の変形ベクトルおよび前記第２の変形ベクトルに基づいて、前記第１の３次元画像上の各点における、前記第１の３次元画像および前記第２の３次元画像間の相対的な磁場歪みを表す磁場歪みベクトルを算出することを特徴とする磁場歪み算出方法。
8. 被検体の対象部位をＭＲＩ装置により撮影にすることにより取得された、撮影時期が異なる第１の３次元画像および第２の３次元画像を取得する手順と、
   前記第１の３次元画像および前記第２の３次元画像のそれぞれから、互いに対応する第１の形状不変領域および第２の形状不変領域を抽出する手順と、
   前記第１の形状不変領域および前記第２の形状不変領域を剛体位置合わせして、前記第１の形状不変領域上の各点の前記第２の形状不変領域上の対応する各点に対する第１の変形ベクトルを取得する手順と、
   前記第１の３次元画像と前記第２の３次元画像とを、前記第１の形状不変領域および前記第２の形状不変領域において非剛体位置合わせして、前記第１の３次元画像の前記第１の形状不変領域上の各点の、前記第２の３次元画像の前記第２の形状不変領域上の対応する各点に対する第２の変形ベクトルを取得する手順と、
   前記第１の変形ベクトルおよび前記第２の変形ベクトルに基づいて、前記第１の３次元画像上の各点における、前記第１の３次元画像および前記第２の３次元画像間の相対的な磁場歪みを表す磁場歪みベクトルを算出する手順とをコンピュータに実行させることを特徴とする磁場歪み算出プログラム。