(57)【要約】【課題】ＭＲタギング画像から対象物の変位（動き）を正確に測定すること【解決手段】タギング画像をＭＲ撮像装置から入力する（Ｓ１１０）。周波数領域へフーリエ変換により画像を変換する（Ｓ１２０）。ピーク値を、対角でない２つ取ってきて（Ｓ１３０）、その逆変換から縞画像を生成する（Ｓ１４０）。この縞画像を基準縞画像として、位相をずらした縞画像を生成する（Ｓ１５０）。  移動後の画像と基準縞画像から得た各位相をずらした縞画像と積をとり、モアレ干渉を生じさせる（Ｓ１６０）。干渉した画像に、縞の周波数より高い周波数部分をカットするために、ローパスフィルタをかける（Ｓ１７０）。このローパスフィルタをかけた画像から、位相画像と、振幅強度画像とを得る（Ｓ１８０）。変位を把握しやすいように、変位を図示することができる（Ｓ２００）。【選択図】図１
2004242796.tif 000002
【特許請求の範囲】
【請求項１】
ＭＲタギング画像による変位測定システムであって、
ＭＲタギング画像と、基準縞画像及び基準縞画像と位相が異なる複数の縞画像との干渉画像を生成する干渉画像生成手段と、
該干渉画像から、位相画像を生成する位相画像生成手段と、
該位相画像から、変位を求める変位計算手段と
を備えることを特徴とする変位測定システム。
【請求項２】
請求項１記載の変位測定システムにおいて、
前記干渉画像生成手段は、ＭＲタギング画像と縞画像との積をとり、ローパスフィルタをかけることで干渉画像を生成することを特徴とする変位測定システム。
【請求項３】
請求項１又は２に記載の変位測定システムにおいて、
さらに、前記干渉画像から、振幅画像を生成する振幅画像生成手段を備えることを特徴とする変位測定システム。
【請求項４】
請求項３に記載の変位測定システムにおいて、
さらに、前記振幅画像から、信頼性の高い変位部分を取り出すためのマスク画像を生成するマスク生成手段を備え、
前記変位計算手段は、前記位相画像をマスク画像で処理してから変位を求めることを特徴とする変位測定システム。
【請求項５】
請求項１〜４のいずれかに記載の変位測定システムにおいて、
前記基準縞画像は、タグ方向ごとにあり、
前記位相の異なる縞画像は、タグ方向ごとにあり、
前記位相画像生成手段及び／又は振幅画像生成手段は、タグ方向ごとに位相画像及び／又は振幅画像を生成し、
前記変位計算手段は、２次元又は３次元の変位を計算する
ことを特徴とする変位測定システム。
【請求項６】
ＭＲタギング画像による変位測定方法であって、
ＭＲタギング画像と、基準縞画像及び基準縞画像と位相が異なる複数の縞画像との干渉画像を生成する干渉画像生成ステップと、
該干渉画像から、位相画像を生成する位相画像生成ステップと、
該位相画像から、変位を求める変位計算ステップと
を備えることを特徴とする変位測定方法。
【請求項７】
請求項６記載の変位測定方法において、
前記干渉画像生成ステップは、ＭＲタギング画像と縞画像との積をとり、ローパスフィルタをかけることで干渉画像を生成することを特徴とする変位測定方法。
【請求項８】
請求項６又は７に記載の変位測定方法において、
さらに、前記干渉画像から、振幅画像を生成する振幅画像生成ステップを備えることを特徴とする変位測定方法。
【請求項９】
請求項８に記載の変位測定方法において、
さらに、前記振幅画像から、信頼性の高い変位部分を取り出すためのマスク画像を生成するマスク生成ステップを備え、
前記変位計算ステップは、前記位相画像をマスク画像で処理してから変位を求めることを特徴とする変位測定方法。
【請求項１０】
請求項６〜９のいずれかに記載の変位測定方法において、
前記基準縞画像は、タグ方向ごとにあり、
前記位相の異なる縞画像は、タグ方向ごとにあり、
前記位相画像生成ステップ及び／又は振幅画像生成ステップは、タグ方向ごとに位相画像及び／又は振幅画像を生成し、
前記変位計算ステップは、２次元又は３次元の変位を計算する
ことを特徴とする変位測定方法。
【請求項１１】
請求項６〜１０のいずれかに記載の変位測定方法をコンピュータ・システムに実行させるコンピュータ・プログラム。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、ＭＲタギング画像から対象物の変位（動き）を測定することに関する。
【０００２】
【技術的背景】
ＭＲ画像は、体内中の原子のスピンを磁場によってそろえ、そのスピンに対してＲＦパルスを照射することで、スピンが放射する共鳴エネルギを画像化したものである。ＭＲタギング画像は、組織に輝度の標識付けを行なったＭＲ画像であり、タギング画像の生成には、ＳＰＡＭＭ  （ｓｐａｔｉａｌ  ｍｏｄｕｌａｔｉｏｎ  ｏｆ  ｍａｇｅｔｉｚａｔｉｏｎ）  が代表的であり、その方法としてＭＲ画像生成の前処理として磁場を空間的に変調し、部分的なスピンの励起状態を生成することで、ＭＲ画像に特定の輝度パターンを付加するものである。
【０００３】
空間的な磁化の変調（ＳＰＡＭＭ）を通じてＭＲＩ中の画像化で二次元的な心筋の動きの画像化に関して画像平面に対して直交する二組のタギング平面を用いることが知られている（特許文献１〜３，非特許文献１，２参照）。
ＭＲＩに関連してタグの特徴を検出する画像処理を用い、配置や歪みに関係する動きマップの中で次に補間（内挿）を用いて特徴を組合せることは知られている（例えば、非特許文献３，４参照のこと）。これらのアプローチはいくつかの手動操作を必要とすることから完全には自動的ではない。
ＭＲタギング画像シーケンスにおいて、オプティカルフローを用いることもまた知られている（非特許文献５〜７参照）。
上記のアプローチでは、箱型関数（ｓｉｎｃ関数を逆フーリエ変換したもの）や、それを代替する飽和タグパターンの代わりに正弦波タグパターンが用いられる。２枚の画像対から得られる時間的な輝度変化と画像の空間的な輝度勾配は、オプティカルフローとして知られる密な動きの推定法を与えるために用いられる特徴である。
そのようなアプローチでは、空間的な輝度の勾配方向に対する一次元の動きに関する情報しか獲得できないことから、動きの算出では正則化手法（近傍の動きを参照すること）が用いられる。
ＭＲタギング画像シーケンスにおいて角度画像を用いることもまた知られている（特許文献４参照）。このアプローチは、一枚の画像から二組のタグの要素を周波数領域で分割して、それぞれのタグの要素に関する角度画像を算出する。連続した二枚の画像でそれぞれの要素に関するオプティカルフローを算出することで、正則化なしに運動を算出することが可能になる。
これらの運動算出に関する問題は、輝度の減衰に伴うタグの消失や、タグが付加していない部分では運動の検出精度が著しく悪いにも関わらず、それらを検出するための方法が存在しないということである。またそれらを手動で検知しなければならないということである。
オプティカルフローの算出では、勾配方程式と呼ばれる対象とする特徴の空間的変化の線形性を必要としており、線形性が十分成り立たないような輝度や角度の空間的な変化の存在する画素での動きの算出では精度が低下する。
同時に時間的な輝度変化が少ない画素では、微小な観測誤差で勾配方程式を満足する領域が広く、算出される動きの信頼性が低下することから、輝度や角度の大きさや振幅の小さな画素での動きの算出では精度が特に低下する。
【０００４】
【特許文献１】
Ｕ．Ｓ．  Ｐａｔ．  Ｎｏ．５０５４４８９
【特許文献２】
Ｕ．Ｓ．  Ｐａｔ．  Ｎｏ．５１１１８２０
【特許文献３】
Ｕ．Ｓ．  Ｐａｔ．  Ｎｏ．５２１７０１６
【特許文献４】
Ｕ．Ｓ．  Ｐａｔ．  Ｎｏ．６４５３１８７
【非特許文献１】
Ａｘｅｌ  ｅｔ  ａｌ．，  ＭＲ  Ｉｍａｇｉｎｇ  ｏｆ  Ｍｏｔｉｏｎ  ｗｉｔｈ  Ｓｐａｔｉａｌ  Ｍｏｄｕｌａｔｉｏｎ  ｏｆ  Ｍａｇｎｅｔｉｚａｔｉｏｎ，  Ｒａｄｉｏｌｏｇｙ，  ｖｏｌ．  １７１，ｐｐ．  ８４１−８４５，  １９８９
【非特許文献２】
Ａｘｅｌ  ｅｔ  ａｌ．，  Ｈｅａｒｔ  Ｗａｌｌ  Ｍｏｔｉｏｎ：  ｉｍｐｒｏｖｅｄ  Ｍｅｔｈｏｄ  ｏｆ  Ｓｐａｔｉａｌ  Ｍｏｄｕｌａｔｉｏｎ  ｏｆＭａｇｎｅｔｉｚａｔｉｏｎ  ｆｏｒ  ＭＲ  ｉｍａｇｉｎｇ，  Ｒａｄｉｏｌｏｇｙ，  ｖｏｌ．  １７２（１），  ３４９−３５０，  １９８９
【非特許文献３】
Ｙｏｕｎｇ  ｅｔ  ａｌ．，  Ｔｈｒｅｅ−ｄｉｍｅｎｓｉｏｎａｌ  ｍｏｔｉｏｎ  ａｎｄ  ｄｅｆｏｒｍａｔｉｏｎ  ｗｉｔｈ  ｓｐａｔｉａｌ  ｍｏｔｉｏｎ  ｏｆ  ｍａｇｎｅｔｉｚａｔｉｏｎ，  Ｒａｄｉｏｌｏｇｙ，  ｖｏｌ．  １８５，  ｐｐ．  ２４１−２４７，  １９９２
【非特許文献４】
ＭｃＶｅｉｇｈ  ｅｔ  ａｌ．，  Ｎｏｎｉｎｖａｓｉｖｅ  ｍｅａｓｕｒｅｍｅｎｔｓ  ｏｆ  ｔｒａｎｓｍｕｒａｌ  ｇｒａｄｉｅｎｔ  ｉｎ  ｍｙｏｃａｒｄｉａｌ  ｓｔｒａｉｎ  ｗｉｔｈ  ＭＲ  ｉｍａｇｉｎｇ，  Ｒａｄｉｏｌｏｇｙ，  ｖｏｌ．  １８０，  ｎｏ．  ３，  ｐｐ．６７７−６８３，  １９９１
【非特許文献５】
Ｐｒｉｎｃｅ  ｅｔ  ａｌ．，  Ｍｏｔｉｏｎ  ｅｓｔｉｍａｔｉｏｎ  ｆｒｏｍ  ｔａｇｇｅｄ  ＭＲ  ｉｍａｇｅ  ｓｅｑｕｅｎｃｅｓ，  ＩＥＥＥ  Ｔｒａｎｓ．  ｏｎ  Ｍｅｄｉｃａｌ  Ｉｍａｇｉｎｇ，  ｖｏｌ．  １１，  ｎｏ．  ２，  ｐｐ．２３８−２４９，  １９９２
【非特許文献６】
Ｇｕｐｔａ  ｅｔ  ａｌ．，  Ｏｎ  ｖａｒｉａｂｌｅ  ｂｒｉｇｈｔｎｅｓｓ  ｏｐｔｉｃａｌ  ｆｌｏｗ  ｆｏｒ  ｔａｇｇｅｄ  ＭＲＩ，  Ｔｅｃｈ．  ｒｅｐ．  ９５−１３，  ＪＨＵ／ＥＣＥ，  １９９５
【非特許文献７】
Ｇｕｐｔａ  ｅｔ  ａｌ．，  Ｂａｎｄｐａｓｓ  ｏｐｔｉｃａｌ  ｆｌｏｗ  ｆｏｒ  ｔａｇｇｅｄ  ＭＲ  Ｉｍａｇｉｎｇ，  Ｐｒｏｃ．  ＩＥＥＥ  ＩＣＩＰ，  １９９７
【０００５】
【発明が解決しようとする課題】
上述の同じ輝度値を探索するような方法では、スピンは時間とともに平衡状態に到達するので、それを画像化したＭＲ画像自体も連続撮像によって輝度が減衰する。同時にタグ付けによってスピンが一方向を向いていない（一様に分布している）部分もまた、磁場の影響を受けることでＭＲ画像本来の輝度に近くなる。これらの作用は、タグの振幅を低下させるように機能することから、どの程度の振幅を持つかということを定量的に評価するための方法が必要となる。また、輝度を追跡する場合には、輝度の変化を考慮した方法が必要になる。
本発明の目的は、ＭＲタギング画像を画像処理することにより、タグから動きを高精度に出すことである。また、タグの存在や残存の強度を数値化し、その動きの信頼度として用いることである。
また、動きを正確に視覚化することも本発明の目的である。
【０００６】
【課題を解決するための手段】
上記目的を達成するために、本発明は、ＭＲタギング画像による変位測定のためのシステムであって、ＭＲタギング画像と、基準縞画像及び基準縞画像と位相が異なる複数の縞画像との干渉画像を生成する干渉画像生成手段と、該干渉画像から、位相画像を生成する位相画像生成手段と、該位相画像から、変位を求める変位計算手段とを備えることを特徴とする。
前記干渉画像生成手段は、ＭＲタギング画像と縞画像との積をとり、ローパスフィルタをかけることで干渉画像を生成することができる。
さらに、前記干渉画像から、振幅画像を生成する振幅画像生成手段を備えることもできる。
さらに、前記振幅画像から、信頼性の高い変位部分を取り出すためのマスク画像を生成するマスク生成手段を備え、前記変位計算手段は、前記位相画像をマスク画像で処理してから変位を求めることができる。
また、前記基準縞画像は、タグ方向ごとにあり、前記位相の異なる縞画像は、タグ方向ごとにあり、前記位相画像生成手段及び／又は振幅画像生成手段は、タグ方向ごとに位相画像及び／又は振幅画像を生成し、前記変位計算手段は、２次元又は３次元の変位を計算することもできる。
その上、上記システムが実行している方法や、上記システムをコンピュータ・システムに構築できるプログラムも本発明である。
【０００７】
【発明の実施の形態】
以下、本発明の実施の形態を、添付図面を参照して説明する。
本発明は、物体の関心領域を空間的に変調するパルスシーケンスを含むＭＲ  画像の動き（変位）を、コンピュータで画像処理することにより、計測するものである。空間的に変調された画像は、タグと同一の基本周波数を有する基準縞と干渉させられることで、空間的に変調した周波数に依存しない位相情報を計算機処理することで獲得する。この位相情報は対象物の変位情報の測定に利用される。同時に位相情報に基づき振幅強度情報を獲得することで、獲得した画素の変位を算出したタグの強度（信頼度）を計算機処理によって定量化する。
【０００８】
図１〜図９を用いて、ＭＲタギング画像（タグ付き画像）から、変位（動き）を求める画像処理について説明する。図１は、本発明の実施形態の画像処理を示すフローチャート、図２は、画像処理過程を示す模式図である。図３から図９は、心臓付近のタグ付き画像の画像処理過程を示す例である。
さて、以下の「位相画像」とは、空間的に変調することでタグ付けされたＭＲ画像に対して干渉を生じさせた際の低周波数成分を抽出した画像の位相を意味し、「振幅強度画像」とはその低周波数成分の個々の画素における振幅の強度を意味する。なお、空間的に変調された画像は、タグの基本周波数と同一の周波数を有する。
【０００９】
まず、タギング画像をＭＲ撮像装置から入力する（Ｓ１１０）。図２（ａ）に示すように、初期状態のタグが乱されていない画像（ｔ＝０）と、変位が欲しい時間でのタギング画像（例えばｔ＝１１など）が必要となる。ここでは、例として、ｔ＝０からの「変位」を算出する。図３は、心臓付近のタグ付き画像の例で、図３（ａ）が初期画像、図３（ｂ）がある時間経過した画像である。
【００１０】
次に、周波数領域へフーリエ変換により画像を変換する（Ｓ１２０）。ｔ＝０の初期画像だけで良い。この周波数領域の画像を図２（ｂ）に示す。ピーク値（ピーク値を含む領域ではない）を、対角でない２つ取ってきて（Ｓ１３０：図２（ｃ）参照）、その逆変換から完全な正弦波画像を生成する（Ｓ１４０：図２（ｄ）参照）。このピークはタグの周波数を表しているので、ピーク値のみ取り出せば、１つの周波数のみ取り出せるので、図２（ｄ）に示すように、単一周期を持つ縞画像が２枚生成できる。これらの画像は、タグの方向にそれぞれ対応している。図４（ａ），図４（ｂ）が図３（ａ）から得た２つの縞画像である。
ピーク値をとるとき、平均分（直流成分）を除けば、１−１ＳＰＡＭＭの二次元画像では４つのピーク（各タグの方向が２つのピークを持つ）が出る。ｓｉｎｃ関数型のＲＦパルスを照射すれば、各象限でタグの基本周波数より高い周波数に幾つかのピークを持つものの、平均分を除けば基本周波数が各象限での最大値を形成するので、この最大値を各象限で探索すれば良い。
この縞画像（図４（ａ），（ｂ））を基準縞画像として、位相をずらした縞画像を生成する（Ｓ１５０）。
【００１１】
移動後の画像（たとえば図３（ｂ））と基準縞画像（図４（ａ）（ｂ））から得た各位相をずらした縞画像と積をとり、モアレ干渉を生じさせる（Ｓ１６０）。この図３（ｂ）と図４（ａ）からの位相をずらした縞画像との積を取った後の画像を図５−１，図５−２に示す。図５−１（ａ）〜図５−２（ｄ）は、それぞれ、位相を１／２πずつずらした縞画像との積をとった画像である。以降は、モアレの原理として後述する手法に従えば良い。
図５−１，図５−２に示した画像に、ローパスフィルタをかけて、縞の周波数より高い周波数部分をカットする（Ｓ１７０）。このときに、干渉縞が形成され、原画像（図３（ａ））とのずれにともない輝度値が変化する画像が得られる。フィルタ後の画像では、縞とのずれが生じていない画素は輝度値が高く、πずれた部分では輝度値が低い。図５−１，図５−２に示した画像に対して、ローパスフィルタをかけた画像を図６−１，図６−２に示す。
このローパスフィルタをかけた画像から、位相画像と、振幅強度画像とを得る（Ｓ１８０）。これらの画像の求め方については、後で詳しく説明する。図６−１（ａ）〜図６−２（ｄ）から求めた位相画像と振幅強度画像とを、それぞれ図７（ａ）と図７（ｂ）に示す。
【００１２】
時間的に推移した場合、タグが部分的に減衰する。または元来タグがほとんど現れない組織も存在することから、これらの評価が必要になる。振幅強度画像では、振幅強度が、元来輝度が小さいか、タグの基本周波数振幅が局所的に小さくなる点では小さくなることから、画像輝度に変位の信頼度としての定量的評価を可能にする。振幅強度画像から、信頼度の高いところのみを抽出するためのマスクを生成する。図７（ｂ）の振幅強度画像から生成したマスク画像を図８（ａ）に、図７（ａ）の位相画像に図８（ａ）のマスクで処理した画像を図８（ｂ）に示す（Ｓ１９０）。
この図８（ｂ）から、変位を把握しやすいように、変位を図示することができる（Ｓ２００）。図８（ｂ）から、心臓部分を切り出したものを図９（ａ）に、この部分の変位を計算して図示したものを図９（ｂ）に示す。図９（ｂ）では、図４に示した２方向の縞画像により求めた変位を合成して、２次元の矢印で示している。
【００１３】
なお、タグの基本周期を既知として始めても良い。このような値は機器に依存する部分なので、最初から知っておけばその周期の画像が撮像可能になるし、いちいち上記の縞画像を求める処理を行なう必要もない。
本発明では、パルスシーケンスとしてＳＰＡＭＭパルスシーケンスを用いることができる。
また、上述では、ＭＲ画像から周波数領域の画像を求めているが、タギング画像を含むＭＲ画像として、ＭＲ装置内にＫスペースと呼ばれる周波数領域での像も格納されている。したがって、それらの装置と直結することによって、画像から周波数領域のスペクトルを抽出する操作を不要とすることが可能である。
すべての画像に対して変位を算出していけば、「動き」を算出することができる。このとき、  ［−π，π）  の周期関数として変位を表現しているため、前から順番に位相を求め、  ［−π，π）    に関する桁上げや桁下げなどの処理が必要となる。
【００１４】
＜詳細なアルゴリズム＞
三次元空間で一般化して（ｉ＝｛１，２，３｝）  説明する。二次元平面の場合、（ｉ＝｛１，２｝）  とすることで適用可能であり、同様に任意の次元について適用可能である。基準縞画像ＴＢｉ，ｉ∈｛１，２，３｝はそれぞれ単一の周波数によって生成される。このとき、基準縞画像上の任意の質点ｐ＝（ｘ，ｙ，ｚ）では、
【数１】
2004242796.tif 000003
と表される。ただし、｛θｉ｝，ｉ∈｛１，２，３｝は初期位相とする。ここで、座標変換
【数２】
2004242796.tif 000004
を行えば、
【数３】
2004242796.tif 000005
と表される。ただし、ＦＦＴが可能なディジタル画像では、｛ｌｉ，ｍｉ，ｎｉ｝，ｉ∈｛１，２，３｝は整数をとる。
【００１５】
タグＴ（ｐ）である
【数４】
2004242796.tif 000006
の個々の要素が１−１ＳＰＡＭＭに代表される単一周波数で生成され、原画像Ｇ０（ｐ）に付加されているとすれば、タグは基準縞画像そのものとなる。すなわち、ｉ∈｛１，２，３｝において、
【数５】
2004242796.tif 000007
である。
【００１６】
一方、タグＴ（ｐ）が箱型関数で生成されている場合、ＤＣ成分を除き基準縞画像の周波数を最小にした整数倍の周波数成分の重ね合わせと見ることができる。すなわち、
【数６】
2004242796.tif 000008
ここで、ｊは有限の周波数成分を表し、｛Ａｉｊ｝は各周波数成分の強度を与える定数であるとする。以降、一般論として１−１ＳＰＡＭＭによるタグも、ｊ＝｛０，１｝とした上式の一つとして議論する。
【００１７】
Ｔ＋１枚の連続ＭＲ画像が撮像されたとし、フレームｔ∈｛０，１，…，Ｔ｝での、タグの付加されていない連続画像の質点ｐｔの輝度をＧｔ（ｐｔ）とする。フレーム０におけるタグ付きの画像Ｉｔ（ｐｔ；｛θｉ｝，｛Ａｉｊ｝）｜ｔ＝０は、Ｇｔ（ｐｔ）｜ｔ＝０に、上記で与えられるタグが付加されたものと見なすことができるから
【数７】
2004242796.tif 000009
と表される。
ｔ≠０の場合、上記関数で与えられるタグが組織の変位・変形にしたがって歪みを持つ。時相０の任意の質点ｐ０が、時相ｔでｐｔ＝ｐ０＋Δｐ，Δｐ＝（Δｕ１，Δｕ２，Δｕ３）Ｔで表される座標に移動するとき、フレームｔにおけるタグのない画像をＧｔ（ｐｔ）とすれば、タギング画像Ｉｔ（ｐｔ）は
【数８】
2004242796.tif 000010
で与えられる。
【００１８】
時相ｔの画像Ｉｔ（ｐｔ；｛θｉ｝，｛Ａｉｊ｝）  と基準縞画像ＴＢｉ（ｕｔｉ；θＢｉ）との積を、Ｍｔｉ（ｐｔ，Δｐｔ，｛θｉ｝，｛Ａｉｊ｝，θＢｉ）とすれば、
【数９】
2004242796.tif 000011
となる。上式において、右辺第１項はバイアス項、第２項はＭＲ画像上に現れたタグそのものを表し、第３項は基準縞画像を表す。また、第４項以降は、個々の周波数成分との積であり、この中で座標ｕｔｉに依存しない項はＡｉ１にかかる第２項だけである。一般にΔｕｔｉはｕｔｉの空間変化よりも小さいから、ＬＰＦによってこの要素のみを抽出する。局所的なＬＰＦ（平均化フィルタ）を適用すれば、Ａｉ０，Ａｉ１も局所的要素として表現できることから、質点ｐｔのそれらをＡｉ０（ｐｔ），Ａｉ１（ｐｔ）として表すことができる。上記の項をＬｔｉとして表し、整理すれば
【数１０】
2004242796.tif 000012
を得る。
【００１９】
位相シフト法は、位相｛θＢｉ｝を変化させた基準縞画像ＴＢｉを用意し、それらの縞との干渉をとることで、干渉縞の位相を得るための方法である。ここで、基準縞画像ＴＢｉの位相θＢｉを２πｋ／Ｋ＋θ１，ｋ∈｛０，…，Ｋ−１｝としたＫ枚の低周波数成分画像Ｌｔｉ（ｐｔ，Δｐｔ；｛θｉ’｝，θＢｉ＝２πｋ／Ｋ）にｓｉｎ（２πｋ／Ｋ）をかけ、ｋに関して足し合わせたものをＳｔｉとすれば、Ｋが４以上の偶数であるとき、
【数１１】
2004242796.tif 000013
となる。同様にｃｏｓ（２πｋ／Ｋ）をかけ、足し合わせたものをＣｔｉとすれば、
【数１２】
2004242796.tif 000014
となる。したがって、Δｕｔｉ（ｐｔ）は、
【数１３】
2004242796.tif 000015
として、求めることができる。ただし、この求めた変位は、［−π，π）の間に限定されている。
【００２０】
時間的に推移した場合、タグが部分的に減衰する。または元来タグがほとんど現れない組織も存在することから、これらの評価が必要になる。このとき、振幅強度Ｗｔｉ（ｐｔ）を
【数１４】
2004242796.tif 000016
として求めれば、元来輝度が小さいか、タグの基本第一周波数振幅Ａｉ１が局所的に小さくなる質点ｐｔではこの値が小さくなることから、画像輝度に変位の信頼度としての定量的評価を可能にする。一般には、タグの存在する方向に関する局所的フィルタリングはできないので、二次元または三次元の局所的ＬＰＦを用いるとすれば、式（９）の注目方向以外のタグも消去され、
【数１５】
2004242796.tif 000017
となる。この場合、同様の処理を行なえば、Δｕｔｉ（ｐｔ），Ｗｔｉ（ｐｔ）は式（１２），（１３）の代わりに、
【数１６】
2004242796.tif 000018
として求められ、注目方向以外のタグの要素が振幅強度に含まれるかどうかの差異が生じる。
マスクを生成する場合、任意のフレームｔに、それぞれの位相に対してＷｔｉ（ｐ）が存在するので、これらを以下のように自乗平均で統合したものを用いるとよい。
【数１７】
2004242796.tif 000019
【００２１】
本方法は、任意の画像（例えば図３（ｂ））のような変位画像から、歪みや変位のない画像（例えば図３（ａ））からみたときの変位の量を直接的に求めることが可能になる。
また、算出される位相の値域が［−π，π）に限定される問題に対して、時相０と時相１の画像を用いて算出した時相１の位相差に対して、時相０と時相２の画像を用いて算出した時相２の位相差が、［−π，π）となるように設定してやることで、任意の時相における位相算出、すなわち変位算出が可能になる。すなわち、同一質点を与えるフレームｔ−１，ｔの座標をｐｔ−１，ｐｔとすれば、式（１２）は
【数１８】
2004242796.tif 000020
と表す。ここで、ｎｔｉ（ｐｔ）はΔｕｔｉ（ｐｔ）の整数部分を表す。この方法では、補間を用いる必要があるが、ｐｔと同じ座標のフレームｔ−１の質点Δｕｔ−１１  （ｐｔ−１）の周期が緩やかであるとすれば、ｐ’ｔ−１によって求められる位相Δｕｔ−１１（ｐ’ｔ−１）によって代替することが可能である。
【００２２】
上述したように、本発明は、小さな変位を求めると同時に、その変位を得るために用いたタグの強度を定量的に求めることで、変位の信頼性を算出するのに用いることができる。
本発明は心臓での使用に特に適しており、虚血性心疾患および心筋梗塞の検出および定量化を評価するのに用いることができる。
また、二つのタグ方向に関する振幅画像から心壁，細胞壁，血管，筋肉，腱，血液，髄液等のいずれかの組合せによる組織の違いを見出すことができる。
使用するＭＲタギング画像としては、グラディエントエコー法，スピンエコー法やスパイラル法からなるグループから選択されたＭＲＩパルスシーケンスを用いることができる。
上述では、図９のように変位を矢印で示した図を示しているが、変位を変位等値線で示すこともできる。
なお、本発明の構成は、上記した実施の形態に限定されるものではなく、本発明の要旨を逸脱しない範囲内において、種々変更を加え得ることは勿論である。
【００２３】
【発明の効果】
以上、説明したように本発明では、ＭＲタギング画像を画像処理することにより、タグから対象物の変位（動き）を高精度に出すことができる。また、タグの存在や残存の強度を数値化し、変位（動き）の信頼度として用いることである。
【図面の簡単な説明】
【図１】本発明の実施形態における処理を示すフローチャートである。
【図２】処理の概略を示す模式図である。
【図３】処理される前の画像例を示す図である。
【図４】２方向の縞画像の例を示す図である。
【図５−１】位相の異なる縞画像との積を取った画像例を示す図である。
【図５−２】位相の異なる縞画像との積を取った画像例の続きを示す図である。
【図６−１】図５−１の画像のローパスフィルタ処理後の画像例を示す図である。
【図６−２】図５−２の画像のローパスフィルタ処理後の画像例を示す図である。
【図７】位相画像と振幅強度画像の例を示す図である。
【図８】振幅強度画像から作成したマスク画像（ａ）と、位相画像に対するマスク処理後の画像（ｂ）を示す図である。
【図９】心臓部分を切り取った位相画像（ａ）と、それから得られた変位を示す画像（ｂ）を示す図である。
【図１】
2004242796.tif 000021
【図２】
2004242796.tif 000022
【図３】
2004242796.tif 000023
【図４】
2004242796.tif 000024
【図５−１】
2004242796.tif 000025
【図５−２】
2004242796.tif 000026
【図６−１】
2004242796.tif 000027
【図６−２】
2004242796.tif 000028
【図７】
2004242796.tif 000029
【図８】
2004242796.tif 000030
【図９】
2004242796.tif 000031