JP4250255B2 - 温度計測方法及び磁気共鳴イメージング装置 - Google Patents
温度計測方法及び磁気共鳴イメージング装置 Download PDFInfo
- Publication number
- JP4250255B2 JP4250255B2 JP11549999A JP11549999A JP4250255B2 JP 4250255 B2 JP4250255 B2 JP 4250255B2 JP 11549999 A JP11549999 A JP 11549999A JP 11549999 A JP11549999 A JP 11549999A JP 4250255 B2 JP4250255 B2 JP 4250255B2
- Authority
- JP
- Japan
- Prior art keywords
- temperature
- subject
- image
- phase
- echo signal
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Images
Landscapes
- Magnetic Resonance Imaging Apparatus (AREA)
Description
【０００１】
【発明の属する技術分野】
本発明は、連続的に被検体中の水素や燐等からの磁気共鳴（以下、MRという）信号を測定し、被検体の温度分布を映像化する磁気共鳴撮影（以下、MRIという）装置を用いた温度計測方法に関する。
【０００２】
【従来の技術】
現在臨床の場で普及しているMRIの撮影対象は、被検体の主たる構成物質であるプロトンである。このプロトンの密度の空間分布や、励起状態の緩和現象の空間分布を画像化することで、人体の頭部、腹部、四肢等の形態又は機能を、２次元もしくは３次元的に撮影する。MRIでの空間分解能は、現状では撮影視野（Field of View；以下、FOVという）当り128、256、512等が用いられている。
【０００３】
近年、ＭＲＩの新しい利用法として、インターベンショナルＭＲＩ（Interventional MRI；以下、IV-MRIと略称する）が普及し始めている。このIV-MRIはオープン型MRI装置の普及とともに広がりつつある。
【０００４】
IV-MRIと他の撮影モダリティによる術中モニタを比較した場合、IV-MRIには、（１）軟部組織の描出能が優れている、（２）X線被曝がなく低侵襲である、（３）任意断面の撮影が可能である、（４）温度モニタが可能である、などの利点がある。
【０００５】
IV-MRIで行われる治療法には、レーザ治療、エタノールなどの薬物注入、高周波（RF）照射切除、超音波治療、低温治療などがある。これらの治療において、MRIの役割は、患部に治療用の穿刺針や細管を挿入する際のリアルタイム・イメージングによるガイド、治療中の組織変化の可視化、加熱・冷却治療中の患部の局所温度のモニタなどである。
【０００６】
IV-MRIに関する参考文献としては下記のものがあげられる。
（１）JF. Schench, FA. Jolesz, PB. Roemer et al; Superconducting open-configuration MR imaging system for imaging-guided therapy; Radiology. Vol. 195, pp805〜814, (1995).
文献（１）には、IV-MRI用のダブルドーナツ型のMRI装置について記載されている。
（２）三井田和夫、原田潤太、土肥美智子、他；インターベンショナルMRIの特徴と問題点；INNERVISION, Vol.12, No.9, pp28〜32, (1997).
文献（２）には、永久磁石型のMRI装置を使用したIV-MRIの基礎的な報告が記載されている。
（３）橋本卓雄、寺尾亨、石橋敏寛、他；MRガイド下経皮的レーザー腰椎椎間板ヘルニア蒸散法；日磁医誌, Vol.18, No.2, pp98〜106, (1998).
文献（３）には、MRモニタを腰椎治療のためのレーザー蒸散術に適用した結果が記載されている。
（４）三井田和夫、原田潤太、土肥美智子、他；オープンタイプMRI装置による透視下のinterventional MRI−特に脳腫瘍生検法について−；日磁医誌, Vol.17, No.8, pp517〜521, (1997).
文献（４）には、永久磁石型のMRI装置を使用したIV-MRIを脳腫瘍生検に適用した結果が記載されている。
【０００７】
IV-MRIにて必要な器具ならびに機能としては、（１）オープン型ガントリ、（２）ガントリ周辺に配置した画像モニタ、（３）専用高周波（RF）受信コイル、（４）リアルタイムイメージング機能、（５）温度変化のモニタリング機能、（６）MR対応の穿刺針などの非磁性器具、などが上げられる。
【０００８】
IV-MRIにて用いられるオープン型MRI装置は、0.2T〜0.5Tの中低磁場の装置で多く実現されている。この理由は、1.5Tなどの高磁場の装置に比べて、中低磁場の装置が、一般にガントリの開放性の点で優れているためである。また、中低磁場でのMRIは、（１）手術器具などによる静磁場の乱れが画像アーチファクトになりにくい、（２）RF照射による生体の加熱が起こりにくい、など、IV-MRIに適した特徴がある。
【０００９】
IV-MRIでは、患部へのアクセス性を向上確保するために、IV-MRI専用のオープン型RF受信コイルが用いられている。このRF受信コイルでは、コイルの導体部分を細くしたり、変形させたりして、穿刺針の患部への挿入が容易になるようにしてある。
【００１０】
IV-MRIでは、目的に応じて撮影シーケンスを使い分けている。すなわち、
（１）病変の広がりについては、通常のMRIシーケンス、
（２）穿刺時のモニタには、GrE（Gradient Echo）系の高速シーケンス、
（３）レーザ治療など、体内の温度をモニタする場合には、信号強度法や位相法（後述する）の温度計測シーケンス、
などが用いられている。また、エタノールなどの薬物治療の場合でも、撮影シーケンスのパラメータを最適化することにより薬物の画像化が可能になる場合が多い。
【００１１】
MRIでは穿刺針による静磁場歪みにより、穿刺針が実際の直径よりも太い陰影として描出される。この理由は、穿刺針による静磁場歪みが穿刺針の周囲にまで及び、信号低下領域が広がるためである。この陰影の広がりの程度は、SE（Spin Echo）シーケンスでは大きい。また、この陰影は、静磁場の方向と穿刺針の方向との関係、撮影シーケンスの読み出し方向によってさまざまに変化する。
【００１２】
MRIでは、生体内の温度をモニタできる。このMRIの機能を活用して、レーザ照射治療のモニタや、RF（Radio Frequency）アブレーションのモニタができるようになりつつある。MRIによる温度モニタの方法は、各種提案されているが、以下に述べる信号強度法と位相法（PPS法：Proton Phase Shift）が多く検討されている。
【００１３】
（１）信号強度法
被検体の温度が変化するとT1値が変化する。また、レーザ蒸散法のように局所的に高温度になる場合、生体のT2値や水分含有量が変化し、信号値が変化する。これらの変化から、生体内の温度変化を定性的に把握することができる。
【００１４】
（２）位相法
プロトンの核磁気共鳴周波数は、温度に対して比例関係にある。共鳴周波数の温度係数Cは、−0.01[ppm/℃]である。GrE系シーケンスで生体の加熱前と加熱後に撮影し、２画像の位相差分を求めると、差分後の位相θ(i,j)[rad]（i,jは画素番号）は、次式で表される。
θ(i,j)＝γ・B0・C・ΔT(i,j)・τ ・・・・・（１）
ここで、γは磁気回転比、B0は静磁場強度、ΔT(i,j)は加熱前後の生体の温度変化、τはエコー時間である。
【００１５】
式（１）を変換することにより、ΔT(i,j)に関する式（２）が得られる。
ΔT(i,j)＝θ(i,j)／(γ・B0・C・τ) ・・・・・（２）
式（２）に基づいて温度変化画像が得られる。温度変化画像の表示には、２次元カラーマップが使われている。
【００１６】
温度画像に関する公知文献としては、例えば、Paul Steiner, Rene Botnar et al; Radio-frequency-induced thermoablation: monitoring with T1-weighted and prpton-frequency-shift MR imaging in an interventional 0.5-T environment; Radiology, Vol.206, pp803〜810, March (1998).がある。この文献では、T1強調法（信号強度法の一種）と位相法を比較し、その優劣について記載している。T1強調法画像、位相法画像はともに２次元画像であり、別途取得した構造画像の上にカラーで温度マップを重ねて表示している。
【００１７】
一方、連続的にMR画像を撮影する際には被検体の体動の影響が問題となるが、この体動除去法に関しては従来ほとんど検討されていなかった。その理由は、従来はMRIの撮影時間が数１０秒と長かったため、連続撮影（Dynamic Scan）はMRI撮影の主流ではなかったためである。
【００１８】
連続撮影の体動除去法の公知例としては、fMRI（functional MRI）などで、オフライン処理によって画像間の位置合わせを行う場合がある。また、一般撮影にて、伸縮性のある幅広の帯を使って、腹部を強く縛り、呼吸動を抑制する方法などがある。
【００１９】
【発明が解決しようとする課題】
生体の加熱時の温度分布を計測するためには、被検体を加熱前と加熱後に撮影し、２画像の位相差分を求めることから、被検体が加熱前後で動くと、温度分布計算結果に誤差が発生する。加熱前後の時間としては最短でも5秒ないし10秒は必要となるし、また加熱に伴う痛みがある場合などもあるので、その間被検体がわずかに動くことは避けられない。また、撮影部位によっては、呼吸や心拍、腸の蠕動などに起因する動きも発生する。このようなわずかな体動がある場合、画像差分時に、とくに臓器の境界部分にアーチファクトが発生する。
【００２０】
しかしながら、従来開示されている体動抑制法は、例えばオフライン処理であったので、患部の温度計測のように、準リアルタイムで温度計測結果を見ながら手術を行うIV-MRIには適用できなかった。また、伸縮性のある幅広の帯を使う方法では、手術部位を覆ってしまうためIV-MRIには適さなかった。
【００２１】
このため、本発明では、MRI装置を用いて被検体内の温度分布を計測する場合に、被検体の動きがあっても被検体内の温度変化を時系列的に安定に計算できるMRI装置を用いた時系列温度計測方法を提供することを目的とする。
【００２２】
【課題を解決するための手段】
上記目的を達成するため、本発明のMRI装置を用いた温度計測方法は、被検体から計測されたエコー時間の異なる第１のエコー信号と第２のエコー信号とを用いて、被検体の温度変化を計測する方法であって、第１のエコー信号と第２のエコー信号とを用いて、第１の複素画像を作成するステップと、第１の複素画像と基準温度を反映した参照複素画像とから、第２の複素画像を作成するステップと、第２の複素画像から位相分布を作成するステップと、位相分布に基づき、被検体の温度変化分布を取得するステップと、を含み、第１のエコー信号はスピンエコー時間より前に計測され、第２のエコー信号は前記スピンエコー時間の後に計測されたエコー信号であることを特徴とする。
本発明の温度計測方法の好ましい一実施形態としては、（１）被検体の温度分布情報を含み、同一の核磁気励起タイミングにて異なるエコー時間を有する複数のMR画像を取得するマルチエコー撮影を連続的に繰り返す時系列マルチエコー撮影を行うステップと、（２）ステップ（１）で得られるマルチエコー複素MR画像のうち、各時相毎に、第１のエコー時間TE1のMR信号から計算された複素MR画像と、第２のエコー時間TE2のMR信号から計算された複素MR画像とから、各時相の第１の複素演算MR画像を作成するステップと、（３）各時相の第１の複素演算MR画像と、基準温度を反映した参照複素演算MR画像とから、各時相の第２の複素演算MR画像を作成するステップと、（４）各時相の第２の複素演算MR画像のアークタンジェントを計算し、位相分布を作成するステップと、（５）ステップ（４）で求めた各時相の位相分布に対して、３次元もしくは２次元の位相アンラップ処理を行うステップと、（６）各時相の位相アンラップ処理後の位相分布に基づき、各時相の３次元もしくは２次元の被検体の温度分布を計算するステップ、とを具備する。
また、上記目的を達成するため、本発明のMRI装置は、被検体からエコー時間の異なる第１のエコー信号と第２のエコー信号の計測を制御する計測制御手段と、第１のエコー信号と第２のエコー信号とを用いて、被検体の温度変化を取得する演算処理手段と、を備え、演算処理手段は、第１のエコー信号と第２のエコー信号とを用いて、第１の複素画像を作成し、第１の複素画像と基準温度を反映した参照複素画像とから、第２の複素画像を作成し、第２の複素画像から位相分布を作成し、位相分布に基づき、被検体の温度変化分布を取得するものであって、 第１のエコー信号はスピンエコー時間より前に計測され、第２のエコー信号は前記スピンエコー時間の後に計測されたエコー信号であることを特徴とする。
【００２３】
この構成では、各時相の被検体の温度分布を計算するための基となる位相分布情報が各時相におけるマルチエコー撮影によって得られた複素MR画像間の演算に基づいて求められたものであるので、両複素MR画像の撮影の間には被検体が動くことは実質的にはなく、本発明の時系列温度計測方法で計測される温度分布では体動の影響を無視することができる。
【００２４】
本発明のMRI装置を用いた時系列温度計測方法では更に、得られた温度分布についてボリュームレンダリング処理を行い、ボリュームレンダリング処理したものを被検体の温度分布画像として表示するステップを具備する。
この構成では、被検体の温度分布画像をボリュームレンダリング法で表示することができるので、レーザ照射などで得られる加熱部分のような局所的に急激な温度変化がある場合でも、時系列で正確に観察することができる。
【００２５】
本発明のMRI装置を用いた時系列温度計測方法では、温度分布画像をボリュームレンダリング処理の元画像データを温度変化領域の広がりに対応して自動抽出するものである。この構成では、温度分布画像上で、温度変化領域のみが画像データの変化する範囲とみて、その温度変化領域の広がりに対応する領域のみボリュームレンダリング処理を行うことにしている。このため、ボリュームレンダリング処理の対象となる元画像データも、温度変化領域の部分に限定することにし、これを自動的に抽出してボリュームレンダリング処理を行うことにより演算処理時間の短縮化を図っている。
【００２６】
【発明の実施の形態】
本発明の実施例を説明する前に、体動がある場合の位相画像の計算誤差について説明する。
本発明に関連して発明者などが検討を行った結果、被検体に体動がある場合の位相画像の計算誤差は、特に臓器の境界領域に強く発生することが明らかになった。これは高周波（RF）磁場のペネトレーション（penetration：貫通）効果によるものと推測される。従って、この効果の影響を除くことが、計算誤差を低減させるのに有効であることがわかった。
【００２７】
ここで、ペネトレーション効果について、図８を用いて簡単に説明する。図８はペネトレーション効果を示す模式図である。図８では、空間に誘電率ε、導電率σの物質（被検体）802があるとする。MRIでは、理想的には空間的に一様なRF磁場801がRF照射コイルから照射される。しかし、被検体802内には、εとσに応じてRF磁場801をキャンセルする方向に渦電流803が発生する。この渦電流803は、ビオサバールの法則により一義的に決定される。
【００２８】
渦電流803が作る誘導RF磁場は、当初のRF磁場801にベクトル加算され、被検体802内のRF磁場804の強度を減らしたり、向きを変えたりする。RF磁場804の向きが局所的に変わるとMR画像の位相も局所的に変化する。従って、異なるεとσを有する臓器間では、MR画像の位相に位相差が生じる。
【００２９】
異なる位相値を有する臓器が含まれる位相画像同士を複素差分する場合、画像間で位置ずれがない理想状態では、これらの位相値はキャンセルされ、ゼロになる。この結果、位相値は温度変化に起因する項のみとなる。一方、現実には体動があるため、臓器の境界部分で、複素差分後の位相値が非ゼロとなり、温度変化に起因する位相値に体動に起因する位相値が加算される。この結果、位相値から計算する温度分布の計算値に誤差が加わることになる。
【００３０】
以下、本発明の実施例を添付画面に沿って説明する。
図６は、本発明の時系列温度計測方法を実施するのに用いるMRI装置の構成例を示すブロック図である。図６において、MRI装置の計測空間600に挿入された被検体601の周囲には、計測空間600に静磁場を発生する磁石602と、計測空間600に傾斜磁場を発生する傾斜磁場コイル603と、計測空間600にRF磁場を発生するRF照射コイル604と、RF磁場に励起されて被検体601が発生するMR信号を検出するRF受信コイル605が配置されている。また、被検体601はベッド612に載置されて、計測空間600に挿入される。
【００３１】
傾斜磁場コイル603は、x軸、y軸、z軸の３方向の傾斜磁場コイルで構成され、傾斜磁場電源609からの信号に応じてそれぞれの方向の傾斜磁場を発生する。RF照射コイル604は、RF送信部610の信号に応じて、RF磁場を発生する。RF受信コイル605が受信したMR信号は、信号検出部606で検出され、信号処理部607で信号処理され、また計算により画像信号に変換される。画像は、表示部608で表示される。
【００３２】
傾斜磁場電源609、RF送信部610、信号検出部606、信号処理部607、表示部608は制御部611で制御される。これらの制御のうち、傾斜磁場電源609、RF送信部610、信号検出部606の制御のタイムチャートは、一般にパルスシーケンスと呼ばれている。
【００３３】
次に、図６のMRI装置を用いて本発明の時系列温度計測方法を実施する手順について、図１のフローチャートに従って説明する。図１は、本発明の時系列温度計測方法の一実施例のフローチャートを示したものである。図示の如く、本実施例の時系列温度計測方法は６個のステップから構成される。以下、図１のフローチャートの各ステップについて順を追って説明する。
【００３４】
先ず、図１の第１ステップ（101）では、被検体601の温度分布情報を含むMR画像データを収集するために時系列マルチエコー撮影シーケンスを実行する。この撮影シーケンスは、被検体601の温度分布情報を含み、同一の核磁気励起タイミングにて異なるエコー時間を有する複数のMR画像を取得するためのもので、連続的に繰り返される。
【００３５】
図２に、本発明に用いる撮影シーケンスの一実施例を示す。この撮影シーケンスは、温度分布情報を含む２次元MR撮影シーケンスであり、マルチスライス２次元オフセットSE法である。本撮影シーケンスでの撮影パラメータは、例えば繰り返し時間TR＝500ms、第１エコー時間TE1＝30ms(＝TE−τ/2)、第２エコー時間TE2＝60ms(＝TE＋τ/2)、撮影マトリックス＝128(読み出し方向)×64(位相エンコード方向)である。また、スライス厚さ＝5mm、積算回数＝１回とする。このとき、撮影時間は32s(＝500×10-3×64s)である。
【００３６】
撮影は少なくとも加熱前と加熱後の２回行う。加熱は、例えばYAGレーザとレーザファイバを使用して、椎間板内にレーザを照射する。１回の照射時間は10秒間とし、所定の治療効果が得られるまでレーザ照射を繰り返す。また、患部の温度変化の時間推移を見るために、２回目の撮影の後、32s毎に撮影を繰り返す。
【００３７】
図２において、本撮影シーケンスでは、90度RFパルス201を印加してから、スピンエコー時間TEの半分TE/2(205)を経過した後に、反転180度RFパルス202を印加し、その後連続して２個のエコーパルス203、204を発生させる。第１のエコーパルス203はスピンエコー時間TEよりτ/2(206)早いタイミングで、第２のエコーパルス204はTEよりτ/2(206)遅いタイミングで発生するように、読み出し傾斜磁場207、208が印加される。図中、209はスライス傾斜磁場、210は位相エンコード傾斜磁場である。
【００３８】
上記の第１のエコー信号203と第２のエコー信号204は、被検体の温度分布情報を含む、MR信号であり、それぞれMR画像に再構成される。ここで得られる２個のMR画像は、第１のエコー時間TE1(≡TE−τ/2)で得られた第１のエコー信号203から計算した複素MR画像Aと、第２のエコー時間TE2(≡TE＋τ/2)で得られた第２のエコー信号204から計算した複素MR画像Bである。
【００３９】
図７に、本発明に用いる撮影シーケンスの他の実施例を示す。この撮影シーケンスは温度分布情報を含む３次元MR撮影シーケンスであり、マルチエコー３次元GrE法である。本撮影シーケンスでの撮影パラメータとしては、例えば、TR＝70ms、TE1＝15ms、TE2＝50ms、撮影マトリックス＝128（読み出し方向）×64（位相エンコード方向）×8（スライスエンコード方向）である。また、バルク厚さ＝32mm、積算回数＝1回とする。このとき、撮影時間は35.8s（＝70×10-3×64×8s）である。
【００４０】
撮影は、少なくとも加熱前と加熱後の２回行う。加熱は、例えばYAGレーザとレーザファイバを使用して、椎間板内にレーザを照射する。１回の照射時間は10秒間とし、所定の治療効果が得られるまでレーザ照射を繰り返す。また、温度の時間的推移を見るためには、２回目の撮影の後、36s毎に撮影を繰り返す。また、上記と同様の撮影は、２次元（2D）マルチスライス撮影でも行うことができる。
【００４１】
図７において、本撮影シーケンスでは、90度RFパルス701を印加後に、読み出し傾斜磁場パルス705、707により、連続して２個のエコーパルス706、708を発生させる。第１のエコーパルス706は、第１のエコー時間TE1(709)のタイミングで、第２のエコーパルス708は第２のエコー時間TE2(710)のタイミングで発生している。上記の第１のエコー信号706と第２のエコー信号708は、被検体の温度分布情報を含むMR信号であり、それぞれMR画像に再構成される。ここで得られる２個のMR画像は、第１のエコー時間TE1で得られた第１のエコー信号706から計算した複素MR画像Aと、第２のエコー時間TE2で得られた第２のエコー信号708から計算した複素MR画像Bである。図中、702はスライス傾斜磁場、703はスライスエンコード傾斜磁場、704は位相エンコード傾斜磁場である。
【００４２】
図１において、第２ステップ（102）から第４ステップ（104）で、第１ステップ（101）で得られた複素MR画像から時系列位相分布を作成する。各ステップの説明の前に、本発明における３次元位相分布の計算について詳述する。
先ず、３次元もしくはマルチスライス２次元PPS法の温度分布画像においては、温度差ΔT(i, j, k)[℃]に対応する位相差θ(i, j, k)[rad]は、式（１）の２次元座標(i, j)を３次元座標(i, j, k)に拡張した式で与えられる。
【００４３】
次に、第１ステップ（101）の撮影で得られた各エコー信号（n：n＝1,2、n＝1の場合は上記の第１のエコー信号に対応し、n＝2の場合は上記の第２のエコー信号に対応する。）を３次元フーリェ変換、もしくはスライス毎に２次元フーリェ変換して、３次元複素MR画像s(i, j, k, n)を得る。ここで、s(i, j, k, n)の絶対値sa(i, j, k, n)＝|s(i, j, k, n)|を求め、所定の閾値以上の画素のみを抽出することによって、被検体が存在する領域を抽出しておく。上記の閾値としては、例えばsa(i, j, k, n)の最大値の20％値を使用する。被検体抽出には、n＝1または2のどちらか一方のエコー信号を用いればよい。
【００４４】
３次元複素MR画像から３次元位相分布を抽出する。先ず、上記の複素MR画像s(i, j, k, n)は、一般に式（３）の如く表される。
s(i, j, k, n)＝sa(i, j, k, n)exp[iφ(i, j, k, n)] ・・・・・（３）
ここで、exp[iφ(i, j, k, n)]は位相成分である。
【００４５】
以下では、上記MR画像は２次元マルチスライスオフセットSEシーケンス法で取得されたものとして説明する。
φ(i, j, k, n)は、一般に、式（４）で記述することができる（特開平5-253192号公報の式（９）など参照）。
φ(i, j, k, n)＝ψ(i, j, k)＋[ΔωB0(i, j, k)＋ΔωT(i, j, k)]×(−1)n・τ/2 ・・（４）
ここで、ψ(i, j, k)はRFペネトレーションによるRF照射磁場位相の空間分布、ΔωB0(i, j, k)は静磁場(B0)の空間的不均一性に伴う核磁気共鳴周波数の変化分、ΔωT(i, j, k)は被検体の温度分布に伴う核磁気共鳴周波数の変化分を示す。
【００４６】
さて、同一撮影シーケンス内の、n＝1とn＝2のMR画像は、同一のRF照射パルスで核を励起し、検出するエコー信号間隔もτ(＝TE2−TE1＝60ms−30ms−30ms)の差である。人体の心拍周期が約1sであることを考慮すると、このτの間に被検体（人体）が動くことは実質的には無く、両画像間で演算する場合は、体動の影響を無視することができる。
【００４７】
そこで、n＝1の画像とn＝2の画像の比ΔSn1/n2(i, j, k)（以下、第１の複素演算画像という）を求めると、
ΔSn1/n2(i, j, k)＝s(i, j, k, 1)／s(i, j, k, 2)
＝[Sa(i, j, k, 1)／Sa(i, j, k, 2)]・exp[i[φ(i, j, k, 1)−φ(i, j, k, 2)]] ・・・（５）
となる。このうちの位相成分exp[i[φ(i, j, k, 1)−φ(i, j, k, 2)]]が温度変化を反映している。式（４）を考慮して、この位相成分を整理すると、
φ(i, j, k, 1)−φ(i, j, k, 2)＝−[ΔωB0(i, j, k)＋ΔωT(i, j, k)]・τ ・・・・・（６）
が得られ、ψ(i, j, k)は除去される。
【００４８】
次に、被検体の加熱前後の複素MR画像の演算について考える。加熱前後の複素演算画像をそれぞれS1(i, j, k)、S2(i, j, k)とし、式（７）、（８）の如く書き直す。
【００４９】
両画像の比ΔS1/2(i, j, k)を求めると。
・（１０）
が得られる。式（１０）においては、指数関数の係数からΔωB0(i, j, k)が除去されている。従って、式（１０）の位相成分
−[ΔωT1(i, j, k)−ΔωT2(i, j, k)]・τ ・・・・・（１１）
は温度変化のみを反映する。
【００５０】
ここで、被検体の加熱前後の局所温度分布をT1(i, j, k)およびT2(i, j, k)とすると、これらの局所温度分布は上記のΔωT1(i, j, k)およびΔωT2(i, j, k)とそれぞれ式（１２）、（１３）の如く関係付けられる。
ΔωT1(i, j, k)＝γ・B0・C・T1(i, j, k) ・・・・・（１２）
ΔωT2(i, j, k)＝γ・B0・C・T2(i, j, k) ・・・・・（１３）
ここで、ΔT(i, j, k)≡T2(i, j, k)−T1(i, j, k)とおくことにより、式（１１）は式（１）を３次元的に記述した下記の式（１a）、すなわちθ(i, j, k)と等価になる。
θ(i, j, k)＝γ・B0・C・ΔT(i, j, k)・τ ・・・・・（１a）
【００５１】
式（９）の演算に使用するデータに関しては、加熱前のものと加熱後のものとでは互いに取得時間が約30秒以上ずれている。従って、その間には被検体が動いている場合が多い。しかしながら、IV-MRIで行われる0.2Tから0.5Tの静磁場強度のオープン型MRI装置を使用した場合、静磁場歪みは静磁場磁石の形状によって決定される場合が殆どであるため、時間変動は少ないと考えてよい。従って、計算結果に誤差を取り込む可能性は小さい。
【００５２】
次に、式（９）で計算されるΔS1/2(i, j, k)の実部と虚部をそれぞれ、
ΔSr1/2(i, j, k)＝real[ΔS1/2(i, j, k)] ・・・・・（１４）
ΔSi1/2(i, j, k)＝img[ΔS1/2(i, j, k)] ・・・・・（１５）
で表すと、その位相成分Δφ1/2(i, j, k)は次式で計算できる。
Δφ1/2(i, j, k)＝arctan[ΔSi1/2(i, j, k)／ΔSr1/2(i, j, k)] ・・・・・（１６）
【００５３】
上記の本発明における３次元位相分布の計算に基づき、第２ステップ（102）から第４ステップ（104）の内容について説明する。
先ず、第２ステップ（102）においては、第１ステップ（101）で得られた第１のエコー信号(203,706)から計算した複素MR画像と第２のエコー信号(204,708)から計算した複素MR画像とから第１の複素演算MR画像を作成する。具体的には、第１ステップ（101）で得られた複素MR画像は式(3)で表されるもの、すなわち、S(i, j, k, 1)が第１のエコー信号から計算した複素MR画像、S(i, j, k, 2)が第２のエコー信号から計算した複素MR画像であり、第１の複素演算MR画像は式（５）で表されるものである。
【００５４】
次に、第３ステップ（103）では、各時相の第１の複素演算MR画像と、基準温度（加熱前温度）を反映した参照複素MR画像とから第２の複素演算MR画像を時系列に作成する。具体的には、加熱前の複素演算MR画像である式（７）のS1(i, j, k)と加熱後の複素演算画像である式（８）のS2(i, j, k)とから、式（９）または（１０）で表される第２の複素演算MR画像であるΔS1/2(i, j, k)を求めるものである。この具体的な例では、加熱前の第１の複素演算MR画像を基準にして、加熱後の第１の複素演算MR画像を時系列的に一定時間間隔（例えば、32秒間隔）で複数枚作成して、それらの各々の第１の複素演算MR画像と加熱前の第１の複素演算MR画像とから各々の第２の複素演算MR画像を作成することにより、加熱後の温度変化に対応する複数枚の時系列的な第２の複素演算MR画像が得られる。ここで、第２の複素演算MR画像を表す式（１０）からは、その位相成分を表す式（１１）が求まり、この式（１１）は加熱前後の温度変化に対応している（式（１２）、（１３）、（１a）参照）。
【００５５】
次に、第４ステップ（104）では、第２の複素演算MR画像のアークタンジェントを計算して、時系列位相分布を作成する。具体的には、式（１４）、（１５）、（１６）に基づいて、第３ステップ（103）で求めた時系列的な第２の複素演算MR画像から、各時相の位相成分Δφ1/2(i, j, k)を計算するものである。
【００５６】
ここで、逆三角関数arctanは、どのような値に対しても−πから＋πまでの値を与えるため、得られた位相差Δφ1-2(i, j, k)の画像はエリアシングを発生する。図３は、これを示したもので、図では見やすくするために２次元表示してある。以上の理由により、Δφ1-2(i, j, k)は、式（１）のθ(i, j, k)には直接対応していない。これを解決するために、次に述べる位相アンラップ処理を行う。
【００５７】
図１のフローチャートにおいて、次の第５ステップ（105）では、第４ステップ（104）で得られた３次元位相分布の計算結果について、３次元位相アンラップ処理を行う。ここで提案する位相アンラップ処理は、レーザ等によって被検体が局所的に加熱された場合、被検体内部の温度分布は連続的であるとの仮定に基づく。この仮定のもとに、３次元画像の被検体部分の非連続点を検出し、非連続点についてはエアリシングが発生していると判断して、位相値のアンラップ処理を行うものである。
【００５８】
以下、位相値のアンラップ処理について図３を参照しながら説明する。図３は本発明に用いる位相アンラップ処理を説明するための図である。位相アンラップ処理は、（１）基点の決定、（２）処理の走査方向の決定、（３）アンラップ処理の実行の手順で行われる。
【００５９】
(1)基点(i0, j0, k0)を決める。
図３の上図は被検体の位相図301を示したものである。この位相図301は２次元のもの（Δφ1-2(i, j, ko)）で、画素i、jについて表示されている。位相図301上には、加熱部分302、基点303、エリアシング304、305が表示されている。先ず、基点(i0, j0, k0)の決定方法としては、被検体内に存在し、かつ温度変化をしていない部位の１画素を選択する。温度変化をしている部位をおおまかに知るためには、温度上昇している領域の信号が低下することを利用する。すなわち、絶対値画像（式（７）、（８）参照）の加熱前後の信号差、
ΔSa1-2(i, j, k, 1)＝| Sa1(i, j, k, 1)−Sa2(i, j, k, 1)|・・・・・（１７）が、最大となる画素を温度変化があった部位と考え、基点から除外する。
【００６０】
基点303の決定については、元になる画像信号値Sa1(i, j, k, 1)およびSa2(i, j, k, 1)が大きいほうが、位相値の精度も高く、基点303として適している。すなわち、被検体内で温度変化（差分画像の絶対値の変化）が少なく、かつ元画像の信号量が大きい画素を基点303とするのが望ましい。
【００６１】
温度変化をしている部位をおおまかに知るための他の方法について図４により説明する。図４も図３の上図と同様に被検体の位相図401を示したものである。位相図401内には、温度変化した領域402、基点403が含まれている。この方法では、図４に示す如く撮影領域を数個（図４では９個）のブロックに分割し、各々のブロックで、被検体内のΔSa1-2(i, j, k, 1)の平均値Smeann（ｎ＝1〜ブロック数）を求め、Smeannが最大となるブロックn内を基点403の選定から除外する。この方法によれば、輝点ノイズなどに起因する不可避的な誤検出を避けることができる。
【００６２】
(2)処理の走査方法を決定する。
アンラップ処理の走査方向は最終的には全方向について行う。しかし、処理の順番は任意でよい。ここでは、先ず基点(i0, j0, k0)303を基準にしてiの正方向を走査方向とする。
【００６３】
(3)アンラップ処理を実行する。
基点(i0, j0, k0)303から出発し、走査方向に沿って順次、隣接画素間の位相差の差分を求め、位相値の差分が位相閾値thを越えていれば、アンラップ処理を行う。位相閾値thとしては、例えばπとする。具体的には、式（１８）に示すような演算を行う。
【００６４】
式(18)の演算内容を図３の下図を用いて説明する。先ず、図３の上図の被検体の位相図301のj＝j0線306上の位相差分Δφ1-2(i, j, k)が図３の下図の鋸歯状の実線Aである。この実線Aは、j＝j0線306上の点Cと点Dに対応する位置で位相値のジャンプをしている。
【００６５】
以下、順にiの正方向にi0+1、i0+2、i0+3、…について、アンラップ処理を繰り返す。被検体の点Eまでこの処理を行ったら、次に基点(i0, j0, k0)303を始点として、iの負方向にi0−1、i0−2、i0−3、…について、θ(i0, j0, k0)を初期値として同様の処理を行う。この結果、画素(i0, j0, k0)の線306についての、Δφ1-2(i, j0, k0)からθ(i, j0, k0)への変換が完了する。
【００６６】
ここで、位相アンラップ処理を正確に行うための前提条件として、「被検体の形状は、凸面で囲まれており、凹部はなく、従って、上記処理を行うにあたり、被検体領域がiの正または負方向に飛び飛びになることはない」がある。これは、通常の臨床利用ではこの条件が満たされていると考えられることから、実用上大きな制約となることはない。
【００６７】
次に、画素座標(i, j0, k0)を初期値として、処理の走査方向を、jもしくはkの方向とし、同様のアンラップ処理を行う。この結果、被検体内部について、３次元的な位相アンラップ処理が完了する。
【００６８】
位相アンラップ処理を更に安定にするためには、直交座標の３方向だけでなく、さらに斜方向、例えばi＝jと平行な方向にも同様なアンラップ処理を行うことが望ましい。
【００６９】
被検体外部の位相値は、背景ノイズから計算される位相値であり、ランダムな値である。これらの位相値は温度分布情報を含まないので、被検体外部の位相値θ(i, j, k)は0とする。
以上の位相アンラップ処理を行った結果、位相画像に存在したエリアシングが除去される。
【００７０】
図１のフローチャートにおいて、次の第６ステップ（106）では、第５ステップ（105）でアンラップ処理した３次元位相分布に基づいて３次元温度分布の計算を行う。すなわち、上記で計算したθ(i, j, k)に基づいて、被検体の温度変化の分布を、式（２）を３次元的表記した式（１２）を使って３次元的に求める。
ΔT(i, j, k)=θ(i, j, k)／(γ・B0・C・τ) ・・・・・（１９）
ここで、γ・B0は共鳴角周波数を与える。B0＝0.3Tの場合には、25.6×106・π[rad/s]である。τは撮影シーケンスのTE値（＝35ms）と等しい。
【００７１】
以上の如く、図１のフローチャートに示した手順に従うことにより、被検体の時系列温度分布を計算することができる。次に、この時系列温度分布の計算結果を表示する方法について説明する。
【００７２】
以下の実施例では、図１のフローチャートの第６ステップで求めた３次元温度分布の計算結果に基づき、表示のための３次元温度画像を計算し、これを表示画面に表示する。本発明においては、この表示に３次元投影法を採用する。温度分布の３次元投影法としては、ボリュームレンダリング法や最大投影値法などがあるが、ボリュームレンダリング法を用いるのが良い。何となれば、レーザによる被検体の加熱はレーザ照射部の中心部分の温度が最も高く、周囲に行くに従って温度が低くなるためである。
【００７３】
これに対し、最大投影値法を用いた場合、基本的には、加熱部分を３次元的に可視化できるものの、レーザファイバ先端部付近の閾値よりも高い高温度部が視覚的にわかりにくいと予想される。このため、本実施例では、温度画像をボリュームレンダリング法で可視化する。
【００７４】
ボリュームレンダリング法では、３次元空間に対してレイキャスティングを行う。この方法としては、例えば、front to back法、back to front法、並行投影法などが当業者に公知である。１画素に注目した場合、αの不透明度をもつ画素にCinの光が入って来たとき、画素の色をCとすると、画素を出て行く光Coutは式（２０）で表される。このCoutが撮影時のボクセル値に相当する。
Cout＝Cin(1−α)＋C・α ・・・・・（２０）
【００７５】
αは0で完全な透明、１で完全な不透明である。αは温度ΔT(i, j, k)に対応する。ΔT(i, j, k)＝0に対しては、α＝0とする。ΔT(i, j, k)の最大値は、例えばレーザ照射では300℃程度であるので、ΔT(i, j, k)＝300[℃]に対して、α＝1とする。すなわち、α＝ΔT(i, j, k)/300とする。このように設定することにより、レーザ照射点が不透明になり、その周辺に、楕円球状に広がる高温度領域を３次元的に画像化することができる。
【００７６】
温度計測のMR撮影を２回以上繰り返す場合には、n回目の撮影と１回目の撮影との間で、上記と同様の演算を行う。その結果、n回目の撮影と１回目の撮影の間に生じた温度変化の画像を３次元的に得ることができる。
【００７７】
さて、ボリュームレンダリング法は有用であるが、演算時間が長いため実用上使い勝手が若干悪いことは周知である。温度画像を連続的に、時系列的に、３次元表示したい場合に演算時間が長いことが問題となる。本実施例の撮影時間は32sまたは35.8sであるが、この数値は、例えばGrE型のマルチショットEPIシーケンスを導入することにより、1/4程度の時間、すなわち10s弱にまで短縮することが可能である。従って、３次元温度分布表示についても10s以内に完了することが重要な課題となる。
【００７８】
以下では、上記課題を解決するためのダイナミック３次元温度計測に適したボリュームレンダリングの改良手法について図５を参照しながら説明する。図５は本発明に用いるボリュームレンダリング処理の一実施例の一部を示す図である。ボリュームレンダリング法の演算時間を減らすには、元になる３次元温度画像の画素数を減らすことが本質的である。局所温度治療では、図５に示す如く、注目する温度変化部分501は、ほぼ球状に局在している。これは、レーザ照射のためのレーザファイバ502やRFトランスジューサが、侵襲的に患部は挿入され、その部位のみを加熱するためである。同様に、焦点超音波による加熱においても、局所のみの加熱となっている。
【００７９】
このように注目点が局在しているということは、血管像を表示するMRAの場合とは大きく異なる。温度画像のこの特徴を利用することによって、元の温度画像からボリュームレンダリングの対象となる画素のみを効率良く自動抽出することができる。
【００８０】
具体的には、温度画像ΔT(i, j, k)が表示される温度表示領域505のうちの、温度変化の最大点もしくは温度変化領域の重心点Q(imax, jmax, kmax)503を求め、この点Q503を中心としてユーザが設定する半径r0の球状領域504もしくは一辺の長さa0の立方体領域（図示せず）のみをボリュームレンダリングの対象領域とする。
【００８１】
球状領域504の半径r0は典型的には20mm程度、立方体領域の一辺の長さa0は40mm程度とする。撮影視野（FOV）が250mm×250mmの画像では、画素サイズは1.95mm(=250mm/128)なので、立方体領域の場合の一辺は20画素(＝40mm/1.95mm)である。また、ボリュームレンダリングに使われる画素数は、20×20×20となる。
【００８２】
ボリュームレンダリング対象領域504の一辺の大きさ（または半径の大きさ）は、撮影毎に変化させても良い。これは、レーザ加熱の場合、最初はごく小さい領域が加熱され、徐々に加熱領域が広がっていくことを考慮して、温度画像ΔT(i, j, k)の表示領域505のうちの、温度変化の最大点Q(imax, jmax, kmax)503を基準としてΔT(i, j, k)の空間的な変化をr(r＝|(i−imax, j−jmax, k−kmax)|)の関数として求める。ΔT(i, j, k)の値が最大変化点の値ΔT(imax, jmax, kmax)の半分になるような半径rを温度変化部位の広がりとして捕らえることができる。これをrdとする。
【００８３】
このときのボリュームレンダリング対象領域504の半径r0の値は、撮影毎に（すなわちボリュームレンダリングの元画像データの組毎に）、例えばr0＝3× rdとなるように変化させる。このように半径r0を設定することにより、ボリュームレンダリング対象の元画像は、温度変化領域の広がりとともに拡大していき、ボリュームレンダリング処理を行うのに必要十分な元画像データから構成され、ボリュームレンダリング処理の演算時間も最短となる。
【００８４】
上記の実施例では、レーザによる加熱時の被検体の温度変化を３次元的に可視化する例を示したが、他の実施例として、例えば、RFアブレーション時の被検体の温度変化を画像化することもできる。また、例えばクライオサージェリーでは、−40℃まで患部を冷却した場合、冷凍領域の直径は20mm程度に達するので、このような低温治療時の温度モニタにも本発明は利用することができる。この場合、ボリュームレンダリングにおいて、温度が下がるほどαの値を大きくするように変換式の符号を反転する必要がある。
【００８５】
【発明の効果】
以上説明した如く、本発明によれば、MRI装置を用いて被検体内の温度分布を計測する場合に、被検体の動きがあっても、被検体内の温度変化を安定に計算することができる。
【００８６】
また、従来生体の加熱時の温度分布を計測するためには、加熱前と加熱後に撮影し、２画像の位相差分を求める必要があり、この間に被検体が動くと計算結果に誤差を発生しやすかった。しかし、本発明によれば、加熱前後の時間が5秒ないし10秒必要となる場合や加熱に伴う痛みがある場合などに、被検体が動くことがあったり、呼吸や心拍、腸の蠕動などに起因する動きが発生しても、上記の画像差分時に臓器の境界部分に生じていたアーチファクトを低減することができる。
【００８７】
また、本発明によれば、従来報告されているオフライン処理の体動抑制法と異なり、オンラインで体動アーチファクトを低減することができる。その結果、温度計測のように、準リアルタイムで結果を見ながら手術を行うIV-MRIに好適な画像が得られる。また、伸縮性の幅広の帯を使う方法と異なり、手術部位を覆ってしまうこともないのでIV-MRIに適している。
【００８８】
また、従来はレーザ照射などで得られる加熱部分を３次元で計測し、３次元で表示する手法は検討されていなかったが、本発明によれば、レーザ照射時の局所的に急激な位相変化がある場合に、立体構造を安定かつ正確に表示することができる。また、位相から血管像を求めるPC（Phase contrast）法のように撮影シーケンスのパラメータを、最適にセットする必要もなく、３次元温度計測が可能である。
【００８９】
また、本発明によれば、加熱部分の時系列的な変化を準リアルタイムで３次元的に可視化できるMRI装置を提供でき、加熱治療時の３次元的な温度変化の様子を直ちにチェックすることができる。
【図面の簡単な説明】
【図１】本発明の時系列温度計測方法の一実施例のフローチャート。
【図２】本発明に用いる撮影シーケンスの一実施例。
【図３】本発明に用いるエリアシング処理を説明するための図。
【図４】本発明に用いる位相アンラップ処理を説明するための他の図。
【図５】本発明に用いるボリュームレンダリング処理の一実施例の一部を示す図。
【図６】本発明を実施するのに用いるMRI装置の構成例を示すブロック図。
【図７】本発明に用いる撮影シーケンスの他の実施例。
【図８】ペネトレーション効果を示す模式図。
【符号の説明】
101…第１ステップ
102…第２ステップ
103…第３ステップ
104…第４ステップ
105…第５ステップ
106…第６ステップ
201, 701…90度RFパルス
202…180度RFパルス
203, 706…第１のエコーパルス
204, 708…第２のエコーパルス
205…TE／2
206…τ／2
207, 208, 705, 707…読み出し傾斜磁場
209, 702…スライス傾斜磁場
210, 704…位相エンコード傾斜磁場
301, 401…位相図
302, 402…加熱部分
303, 403…基点
304, 305…エリアシング
306…j＝j0線
501…温度変化部分
502…レーザファイバ
503…温度変化の最大点もしくは温度変化領域の重心点Q
504…球状領域(ボリユームレンダリング対象領域)
505…温度表示領域
600…計測空間
601…被検体
602…磁石
603…傾斜磁場コイル
604…RF照射コイル
605…RF受信コイル
606…信号検出部
607…信号処理部
608…表示部
609…傾斜磁場電源
610…RF送信部
611…制御部
612…ベッド
703…スライスエンコード傾斜磁場
709…第１のエコー時間
710…第２のエコー時間
801, 804…RF磁場
802…被検体
803…渦電流
【発明の属する技術分野】
本発明は、連続的に被検体中の水素や燐等からの磁気共鳴（以下、MRという）信号を測定し、被検体の温度分布を映像化する磁気共鳴撮影（以下、MRIという）装置を用いた温度計測方法に関する。
【０００２】
【従来の技術】
現在臨床の場で普及しているMRIの撮影対象は、被検体の主たる構成物質であるプロトンである。このプロトンの密度の空間分布や、励起状態の緩和現象の空間分布を画像化することで、人体の頭部、腹部、四肢等の形態又は機能を、２次元もしくは３次元的に撮影する。MRIでの空間分解能は、現状では撮影視野（Field of View；以下、FOVという）当り128、256、512等が用いられている。
【０００３】
近年、ＭＲＩの新しい利用法として、インターベンショナルＭＲＩ（Interventional MRI；以下、IV-MRIと略称する）が普及し始めている。このIV-MRIはオープン型MRI装置の普及とともに広がりつつある。
【０００４】
IV-MRIと他の撮影モダリティによる術中モニタを比較した場合、IV-MRIには、（１）軟部組織の描出能が優れている、（２）X線被曝がなく低侵襲である、（３）任意断面の撮影が可能である、（４）温度モニタが可能である、などの利点がある。
【０００５】
IV-MRIで行われる治療法には、レーザ治療、エタノールなどの薬物注入、高周波（RF）照射切除、超音波治療、低温治療などがある。これらの治療において、MRIの役割は、患部に治療用の穿刺針や細管を挿入する際のリアルタイム・イメージングによるガイド、治療中の組織変化の可視化、加熱・冷却治療中の患部の局所温度のモニタなどである。
【０００６】
IV-MRIに関する参考文献としては下記のものがあげられる。
（１）JF. Schench, FA. Jolesz, PB. Roemer et al; Superconducting open-configuration MR imaging system for imaging-guided therapy; Radiology. Vol. 195, pp805〜814, (1995).
文献（１）には、IV-MRI用のダブルドーナツ型のMRI装置について記載されている。
（２）三井田和夫、原田潤太、土肥美智子、他；インターベンショナルMRIの特徴と問題点；INNERVISION, Vol.12, No.9, pp28〜32, (1997).
文献（２）には、永久磁石型のMRI装置を使用したIV-MRIの基礎的な報告が記載されている。
（３）橋本卓雄、寺尾亨、石橋敏寛、他；MRガイド下経皮的レーザー腰椎椎間板ヘルニア蒸散法；日磁医誌, Vol.18, No.2, pp98〜106, (1998).
文献（３）には、MRモニタを腰椎治療のためのレーザー蒸散術に適用した結果が記載されている。
（４）三井田和夫、原田潤太、土肥美智子、他；オープンタイプMRI装置による透視下のinterventional MRI−特に脳腫瘍生検法について−；日磁医誌, Vol.17, No.8, pp517〜521, (1997).
文献（４）には、永久磁石型のMRI装置を使用したIV-MRIを脳腫瘍生検に適用した結果が記載されている。
【０００７】
IV-MRIにて必要な器具ならびに機能としては、（１）オープン型ガントリ、（２）ガントリ周辺に配置した画像モニタ、（３）専用高周波（RF）受信コイル、（４）リアルタイムイメージング機能、（５）温度変化のモニタリング機能、（６）MR対応の穿刺針などの非磁性器具、などが上げられる。
【０００８】
IV-MRIにて用いられるオープン型MRI装置は、0.2T〜0.5Tの中低磁場の装置で多く実現されている。この理由は、1.5Tなどの高磁場の装置に比べて、中低磁場の装置が、一般にガントリの開放性の点で優れているためである。また、中低磁場でのMRIは、（１）手術器具などによる静磁場の乱れが画像アーチファクトになりにくい、（２）RF照射による生体の加熱が起こりにくい、など、IV-MRIに適した特徴がある。
【０００９】
IV-MRIでは、患部へのアクセス性を向上確保するために、IV-MRI専用のオープン型RF受信コイルが用いられている。このRF受信コイルでは、コイルの導体部分を細くしたり、変形させたりして、穿刺針の患部への挿入が容易になるようにしてある。
【００１０】
IV-MRIでは、目的に応じて撮影シーケンスを使い分けている。すなわち、
（１）病変の広がりについては、通常のMRIシーケンス、
（２）穿刺時のモニタには、GrE（Gradient Echo）系の高速シーケンス、
（３）レーザ治療など、体内の温度をモニタする場合には、信号強度法や位相法（後述する）の温度計測シーケンス、
などが用いられている。また、エタノールなどの薬物治療の場合でも、撮影シーケンスのパラメータを最適化することにより薬物の画像化が可能になる場合が多い。
【００１１】
MRIでは穿刺針による静磁場歪みにより、穿刺針が実際の直径よりも太い陰影として描出される。この理由は、穿刺針による静磁場歪みが穿刺針の周囲にまで及び、信号低下領域が広がるためである。この陰影の広がりの程度は、SE（Spin Echo）シーケンスでは大きい。また、この陰影は、静磁場の方向と穿刺針の方向との関係、撮影シーケンスの読み出し方向によってさまざまに変化する。
【００１２】
MRIでは、生体内の温度をモニタできる。このMRIの機能を活用して、レーザ照射治療のモニタや、RF（Radio Frequency）アブレーションのモニタができるようになりつつある。MRIによる温度モニタの方法は、各種提案されているが、以下に述べる信号強度法と位相法（PPS法：Proton Phase Shift）が多く検討されている。
【００１３】
（１）信号強度法
被検体の温度が変化するとT1値が変化する。また、レーザ蒸散法のように局所的に高温度になる場合、生体のT2値や水分含有量が変化し、信号値が変化する。これらの変化から、生体内の温度変化を定性的に把握することができる。
【００１４】
（２）位相法
プロトンの核磁気共鳴周波数は、温度に対して比例関係にある。共鳴周波数の温度係数Cは、−0.01[ppm/℃]である。GrE系シーケンスで生体の加熱前と加熱後に撮影し、２画像の位相差分を求めると、差分後の位相θ(i,j)[rad]（i,jは画素番号）は、次式で表される。
θ(i,j)＝γ・B0・C・ΔT(i,j)・τ ・・・・・（１）
ここで、γは磁気回転比、B0は静磁場強度、ΔT(i,j)は加熱前後の生体の温度変化、τはエコー時間である。
【００１５】
式（１）を変換することにより、ΔT(i,j)に関する式（２）が得られる。
ΔT(i,j)＝θ(i,j)／(γ・B0・C・τ) ・・・・・（２）
式（２）に基づいて温度変化画像が得られる。温度変化画像の表示には、２次元カラーマップが使われている。
【００１６】
温度画像に関する公知文献としては、例えば、Paul Steiner, Rene Botnar et al; Radio-frequency-induced thermoablation: monitoring with T1-weighted and prpton-frequency-shift MR imaging in an interventional 0.5-T environment; Radiology, Vol.206, pp803〜810, March (1998).がある。この文献では、T1強調法（信号強度法の一種）と位相法を比較し、その優劣について記載している。T1強調法画像、位相法画像はともに２次元画像であり、別途取得した構造画像の上にカラーで温度マップを重ねて表示している。
【００１７】
一方、連続的にMR画像を撮影する際には被検体の体動の影響が問題となるが、この体動除去法に関しては従来ほとんど検討されていなかった。その理由は、従来はMRIの撮影時間が数１０秒と長かったため、連続撮影（Dynamic Scan）はMRI撮影の主流ではなかったためである。
【００１８】
連続撮影の体動除去法の公知例としては、fMRI（functional MRI）などで、オフライン処理によって画像間の位置合わせを行う場合がある。また、一般撮影にて、伸縮性のある幅広の帯を使って、腹部を強く縛り、呼吸動を抑制する方法などがある。
【００１９】
【発明が解決しようとする課題】
生体の加熱時の温度分布を計測するためには、被検体を加熱前と加熱後に撮影し、２画像の位相差分を求めることから、被検体が加熱前後で動くと、温度分布計算結果に誤差が発生する。加熱前後の時間としては最短でも5秒ないし10秒は必要となるし、また加熱に伴う痛みがある場合などもあるので、その間被検体がわずかに動くことは避けられない。また、撮影部位によっては、呼吸や心拍、腸の蠕動などに起因する動きも発生する。このようなわずかな体動がある場合、画像差分時に、とくに臓器の境界部分にアーチファクトが発生する。
【００２０】
しかしながら、従来開示されている体動抑制法は、例えばオフライン処理であったので、患部の温度計測のように、準リアルタイムで温度計測結果を見ながら手術を行うIV-MRIには適用できなかった。また、伸縮性のある幅広の帯を使う方法では、手術部位を覆ってしまうためIV-MRIには適さなかった。
【００２１】
このため、本発明では、MRI装置を用いて被検体内の温度分布を計測する場合に、被検体の動きがあっても被検体内の温度変化を時系列的に安定に計算できるMRI装置を用いた時系列温度計測方法を提供することを目的とする。
【００２２】
【課題を解決するための手段】
上記目的を達成するため、本発明のMRI装置を用いた温度計測方法は、被検体から計測されたエコー時間の異なる第１のエコー信号と第２のエコー信号とを用いて、被検体の温度変化を計測する方法であって、第１のエコー信号と第２のエコー信号とを用いて、第１の複素画像を作成するステップと、第１の複素画像と基準温度を反映した参照複素画像とから、第２の複素画像を作成するステップと、第２の複素画像から位相分布を作成するステップと、位相分布に基づき、被検体の温度変化分布を取得するステップと、を含み、第１のエコー信号はスピンエコー時間より前に計測され、第２のエコー信号は前記スピンエコー時間の後に計測されたエコー信号であることを特徴とする。
本発明の温度計測方法の好ましい一実施形態としては、（１）被検体の温度分布情報を含み、同一の核磁気励起タイミングにて異なるエコー時間を有する複数のMR画像を取得するマルチエコー撮影を連続的に繰り返す時系列マルチエコー撮影を行うステップと、（２）ステップ（１）で得られるマルチエコー複素MR画像のうち、各時相毎に、第１のエコー時間TE1のMR信号から計算された複素MR画像と、第２のエコー時間TE2のMR信号から計算された複素MR画像とから、各時相の第１の複素演算MR画像を作成するステップと、（３）各時相の第１の複素演算MR画像と、基準温度を反映した参照複素演算MR画像とから、各時相の第２の複素演算MR画像を作成するステップと、（４）各時相の第２の複素演算MR画像のアークタンジェントを計算し、位相分布を作成するステップと、（５）ステップ（４）で求めた各時相の位相分布に対して、３次元もしくは２次元の位相アンラップ処理を行うステップと、（６）各時相の位相アンラップ処理後の位相分布に基づき、各時相の３次元もしくは２次元の被検体の温度分布を計算するステップ、とを具備する。
また、上記目的を達成するため、本発明のMRI装置は、被検体からエコー時間の異なる第１のエコー信号と第２のエコー信号の計測を制御する計測制御手段と、第１のエコー信号と第２のエコー信号とを用いて、被検体の温度変化を取得する演算処理手段と、を備え、演算処理手段は、第１のエコー信号と第２のエコー信号とを用いて、第１の複素画像を作成し、第１の複素画像と基準温度を反映した参照複素画像とから、第２の複素画像を作成し、第２の複素画像から位相分布を作成し、位相分布に基づき、被検体の温度変化分布を取得するものであって、 第１のエコー信号はスピンエコー時間より前に計測され、第２のエコー信号は前記スピンエコー時間の後に計測されたエコー信号であることを特徴とする。
【００２３】
この構成では、各時相の被検体の温度分布を計算するための基となる位相分布情報が各時相におけるマルチエコー撮影によって得られた複素MR画像間の演算に基づいて求められたものであるので、両複素MR画像の撮影の間には被検体が動くことは実質的にはなく、本発明の時系列温度計測方法で計測される温度分布では体動の影響を無視することができる。
【００２４】
本発明のMRI装置を用いた時系列温度計測方法では更に、得られた温度分布についてボリュームレンダリング処理を行い、ボリュームレンダリング処理したものを被検体の温度分布画像として表示するステップを具備する。
この構成では、被検体の温度分布画像をボリュームレンダリング法で表示することができるので、レーザ照射などで得られる加熱部分のような局所的に急激な温度変化がある場合でも、時系列で正確に観察することができる。
【００２５】
本発明のMRI装置を用いた時系列温度計測方法では、温度分布画像をボリュームレンダリング処理の元画像データを温度変化領域の広がりに対応して自動抽出するものである。この構成では、温度分布画像上で、温度変化領域のみが画像データの変化する範囲とみて、その温度変化領域の広がりに対応する領域のみボリュームレンダリング処理を行うことにしている。このため、ボリュームレンダリング処理の対象となる元画像データも、温度変化領域の部分に限定することにし、これを自動的に抽出してボリュームレンダリング処理を行うことにより演算処理時間の短縮化を図っている。
【００２６】
【発明の実施の形態】
本発明の実施例を説明する前に、体動がある場合の位相画像の計算誤差について説明する。
本発明に関連して発明者などが検討を行った結果、被検体に体動がある場合の位相画像の計算誤差は、特に臓器の境界領域に強く発生することが明らかになった。これは高周波（RF）磁場のペネトレーション（penetration：貫通）効果によるものと推測される。従って、この効果の影響を除くことが、計算誤差を低減させるのに有効であることがわかった。
【００２７】
ここで、ペネトレーション効果について、図８を用いて簡単に説明する。図８はペネトレーション効果を示す模式図である。図８では、空間に誘電率ε、導電率σの物質（被検体）802があるとする。MRIでは、理想的には空間的に一様なRF磁場801がRF照射コイルから照射される。しかし、被検体802内には、εとσに応じてRF磁場801をキャンセルする方向に渦電流803が発生する。この渦電流803は、ビオサバールの法則により一義的に決定される。
【００２８】
渦電流803が作る誘導RF磁場は、当初のRF磁場801にベクトル加算され、被検体802内のRF磁場804の強度を減らしたり、向きを変えたりする。RF磁場804の向きが局所的に変わるとMR画像の位相も局所的に変化する。従って、異なるεとσを有する臓器間では、MR画像の位相に位相差が生じる。
【００２９】
異なる位相値を有する臓器が含まれる位相画像同士を複素差分する場合、画像間で位置ずれがない理想状態では、これらの位相値はキャンセルされ、ゼロになる。この結果、位相値は温度変化に起因する項のみとなる。一方、現実には体動があるため、臓器の境界部分で、複素差分後の位相値が非ゼロとなり、温度変化に起因する位相値に体動に起因する位相値が加算される。この結果、位相値から計算する温度分布の計算値に誤差が加わることになる。
【００３０】
以下、本発明の実施例を添付画面に沿って説明する。
図６は、本発明の時系列温度計測方法を実施するのに用いるMRI装置の構成例を示すブロック図である。図６において、MRI装置の計測空間600に挿入された被検体601の周囲には、計測空間600に静磁場を発生する磁石602と、計測空間600に傾斜磁場を発生する傾斜磁場コイル603と、計測空間600にRF磁場を発生するRF照射コイル604と、RF磁場に励起されて被検体601が発生するMR信号を検出するRF受信コイル605が配置されている。また、被検体601はベッド612に載置されて、計測空間600に挿入される。
【００３１】
傾斜磁場コイル603は、x軸、y軸、z軸の３方向の傾斜磁場コイルで構成され、傾斜磁場電源609からの信号に応じてそれぞれの方向の傾斜磁場を発生する。RF照射コイル604は、RF送信部610の信号に応じて、RF磁場を発生する。RF受信コイル605が受信したMR信号は、信号検出部606で検出され、信号処理部607で信号処理され、また計算により画像信号に変換される。画像は、表示部608で表示される。
【００３２】
傾斜磁場電源609、RF送信部610、信号検出部606、信号処理部607、表示部608は制御部611で制御される。これらの制御のうち、傾斜磁場電源609、RF送信部610、信号検出部606の制御のタイムチャートは、一般にパルスシーケンスと呼ばれている。
【００３３】
次に、図６のMRI装置を用いて本発明の時系列温度計測方法を実施する手順について、図１のフローチャートに従って説明する。図１は、本発明の時系列温度計測方法の一実施例のフローチャートを示したものである。図示の如く、本実施例の時系列温度計測方法は６個のステップから構成される。以下、図１のフローチャートの各ステップについて順を追って説明する。
【００３４】
先ず、図１の第１ステップ（101）では、被検体601の温度分布情報を含むMR画像データを収集するために時系列マルチエコー撮影シーケンスを実行する。この撮影シーケンスは、被検体601の温度分布情報を含み、同一の核磁気励起タイミングにて異なるエコー時間を有する複数のMR画像を取得するためのもので、連続的に繰り返される。
【００３５】
図２に、本発明に用いる撮影シーケンスの一実施例を示す。この撮影シーケンスは、温度分布情報を含む２次元MR撮影シーケンスであり、マルチスライス２次元オフセットSE法である。本撮影シーケンスでの撮影パラメータは、例えば繰り返し時間TR＝500ms、第１エコー時間TE1＝30ms(＝TE−τ/2)、第２エコー時間TE2＝60ms(＝TE＋τ/2)、撮影マトリックス＝128(読み出し方向)×64(位相エンコード方向)である。また、スライス厚さ＝5mm、積算回数＝１回とする。このとき、撮影時間は32s(＝500×10-3×64s)である。
【００３６】
撮影は少なくとも加熱前と加熱後の２回行う。加熱は、例えばYAGレーザとレーザファイバを使用して、椎間板内にレーザを照射する。１回の照射時間は10秒間とし、所定の治療効果が得られるまでレーザ照射を繰り返す。また、患部の温度変化の時間推移を見るために、２回目の撮影の後、32s毎に撮影を繰り返す。
【００３７】
図２において、本撮影シーケンスでは、90度RFパルス201を印加してから、スピンエコー時間TEの半分TE/2(205)を経過した後に、反転180度RFパルス202を印加し、その後連続して２個のエコーパルス203、204を発生させる。第１のエコーパルス203はスピンエコー時間TEよりτ/2(206)早いタイミングで、第２のエコーパルス204はTEよりτ/2(206)遅いタイミングで発生するように、読み出し傾斜磁場207、208が印加される。図中、209はスライス傾斜磁場、210は位相エンコード傾斜磁場である。
【００３８】
上記の第１のエコー信号203と第２のエコー信号204は、被検体の温度分布情報を含む、MR信号であり、それぞれMR画像に再構成される。ここで得られる２個のMR画像は、第１のエコー時間TE1(≡TE−τ/2)で得られた第１のエコー信号203から計算した複素MR画像Aと、第２のエコー時間TE2(≡TE＋τ/2)で得られた第２のエコー信号204から計算した複素MR画像Bである。
【００３９】
図７に、本発明に用いる撮影シーケンスの他の実施例を示す。この撮影シーケンスは温度分布情報を含む３次元MR撮影シーケンスであり、マルチエコー３次元GrE法である。本撮影シーケンスでの撮影パラメータとしては、例えば、TR＝70ms、TE1＝15ms、TE2＝50ms、撮影マトリックス＝128（読み出し方向）×64（位相エンコード方向）×8（スライスエンコード方向）である。また、バルク厚さ＝32mm、積算回数＝1回とする。このとき、撮影時間は35.8s（＝70×10-3×64×8s）である。
【００４０】
撮影は、少なくとも加熱前と加熱後の２回行う。加熱は、例えばYAGレーザとレーザファイバを使用して、椎間板内にレーザを照射する。１回の照射時間は10秒間とし、所定の治療効果が得られるまでレーザ照射を繰り返す。また、温度の時間的推移を見るためには、２回目の撮影の後、36s毎に撮影を繰り返す。また、上記と同様の撮影は、２次元（2D）マルチスライス撮影でも行うことができる。
【００４１】
図７において、本撮影シーケンスでは、90度RFパルス701を印加後に、読み出し傾斜磁場パルス705、707により、連続して２個のエコーパルス706、708を発生させる。第１のエコーパルス706は、第１のエコー時間TE1(709)のタイミングで、第２のエコーパルス708は第２のエコー時間TE2(710)のタイミングで発生している。上記の第１のエコー信号706と第２のエコー信号708は、被検体の温度分布情報を含むMR信号であり、それぞれMR画像に再構成される。ここで得られる２個のMR画像は、第１のエコー時間TE1で得られた第１のエコー信号706から計算した複素MR画像Aと、第２のエコー時間TE2で得られた第２のエコー信号708から計算した複素MR画像Bである。図中、702はスライス傾斜磁場、703はスライスエンコード傾斜磁場、704は位相エンコード傾斜磁場である。
【００４２】
図１において、第２ステップ（102）から第４ステップ（104）で、第１ステップ（101）で得られた複素MR画像から時系列位相分布を作成する。各ステップの説明の前に、本発明における３次元位相分布の計算について詳述する。
先ず、３次元もしくはマルチスライス２次元PPS法の温度分布画像においては、温度差ΔT(i, j, k)[℃]に対応する位相差θ(i, j, k)[rad]は、式（１）の２次元座標(i, j)を３次元座標(i, j, k)に拡張した式で与えられる。
【００４３】
次に、第１ステップ（101）の撮影で得られた各エコー信号（n：n＝1,2、n＝1の場合は上記の第１のエコー信号に対応し、n＝2の場合は上記の第２のエコー信号に対応する。）を３次元フーリェ変換、もしくはスライス毎に２次元フーリェ変換して、３次元複素MR画像s(i, j, k, n)を得る。ここで、s(i, j, k, n)の絶対値sa(i, j, k, n)＝|s(i, j, k, n)|を求め、所定の閾値以上の画素のみを抽出することによって、被検体が存在する領域を抽出しておく。上記の閾値としては、例えばsa(i, j, k, n)の最大値の20％値を使用する。被検体抽出には、n＝1または2のどちらか一方のエコー信号を用いればよい。
【００４４】
３次元複素MR画像から３次元位相分布を抽出する。先ず、上記の複素MR画像s(i, j, k, n)は、一般に式（３）の如く表される。
s(i, j, k, n)＝sa(i, j, k, n)exp[iφ(i, j, k, n)] ・・・・・（３）
ここで、exp[iφ(i, j, k, n)]は位相成分である。
【００４５】
以下では、上記MR画像は２次元マルチスライスオフセットSEシーケンス法で取得されたものとして説明する。
φ(i, j, k, n)は、一般に、式（４）で記述することができる（特開平5-253192号公報の式（９）など参照）。
φ(i, j, k, n)＝ψ(i, j, k)＋[ΔωB0(i, j, k)＋ΔωT(i, j, k)]×(−1)n・τ/2 ・・（４）
ここで、ψ(i, j, k)はRFペネトレーションによるRF照射磁場位相の空間分布、ΔωB0(i, j, k)は静磁場(B0)の空間的不均一性に伴う核磁気共鳴周波数の変化分、ΔωT(i, j, k)は被検体の温度分布に伴う核磁気共鳴周波数の変化分を示す。
【００４６】
さて、同一撮影シーケンス内の、n＝1とn＝2のMR画像は、同一のRF照射パルスで核を励起し、検出するエコー信号間隔もτ(＝TE2−TE1＝60ms−30ms−30ms)の差である。人体の心拍周期が約1sであることを考慮すると、このτの間に被検体（人体）が動くことは実質的には無く、両画像間で演算する場合は、体動の影響を無視することができる。
【００４７】
そこで、n＝1の画像とn＝2の画像の比ΔSn1/n2(i, j, k)（以下、第１の複素演算画像という）を求めると、
ΔSn1/n2(i, j, k)＝s(i, j, k, 1)／s(i, j, k, 2)
＝[Sa(i, j, k, 1)／Sa(i, j, k, 2)]・exp[i[φ(i, j, k, 1)−φ(i, j, k, 2)]] ・・・（５）
となる。このうちの位相成分exp[i[φ(i, j, k, 1)−φ(i, j, k, 2)]]が温度変化を反映している。式（４）を考慮して、この位相成分を整理すると、
φ(i, j, k, 1)−φ(i, j, k, 2)＝−[ΔωB0(i, j, k)＋ΔωT(i, j, k)]・τ ・・・・・（６）
が得られ、ψ(i, j, k)は除去される。
【００４８】
次に、被検体の加熱前後の複素MR画像の演算について考える。加熱前後の複素演算画像をそれぞれS1(i, j, k)、S2(i, j, k)とし、式（７）、（８）の如く書き直す。
両画像の比ΔS1/2(i, j, k)を求めると。
が得られる。式（１０）においては、指数関数の係数からΔωB0(i, j, k)が除去されている。従って、式（１０）の位相成分
−[ΔωT1(i, j, k)−ΔωT2(i, j, k)]・τ ・・・・・（１１）
は温度変化のみを反映する。
【００５０】
ここで、被検体の加熱前後の局所温度分布をT1(i, j, k)およびT2(i, j, k)とすると、これらの局所温度分布は上記のΔωT1(i, j, k)およびΔωT2(i, j, k)とそれぞれ式（１２）、（１３）の如く関係付けられる。
ΔωT1(i, j, k)＝γ・B0・C・T1(i, j, k) ・・・・・（１２）
ΔωT2(i, j, k)＝γ・B0・C・T2(i, j, k) ・・・・・（１３）
ここで、ΔT(i, j, k)≡T2(i, j, k)−T1(i, j, k)とおくことにより、式（１１）は式（１）を３次元的に記述した下記の式（１a）、すなわちθ(i, j, k)と等価になる。
θ(i, j, k)＝γ・B0・C・ΔT(i, j, k)・τ ・・・・・（１a）
【００５１】
式（９）の演算に使用するデータに関しては、加熱前のものと加熱後のものとでは互いに取得時間が約30秒以上ずれている。従って、その間には被検体が動いている場合が多い。しかしながら、IV-MRIで行われる0.2Tから0.5Tの静磁場強度のオープン型MRI装置を使用した場合、静磁場歪みは静磁場磁石の形状によって決定される場合が殆どであるため、時間変動は少ないと考えてよい。従って、計算結果に誤差を取り込む可能性は小さい。
【００５２】
次に、式（９）で計算されるΔS1/2(i, j, k)の実部と虚部をそれぞれ、
ΔSr1/2(i, j, k)＝real[ΔS1/2(i, j, k)] ・・・・・（１４）
ΔSi1/2(i, j, k)＝img[ΔS1/2(i, j, k)] ・・・・・（１５）
で表すと、その位相成分Δφ1/2(i, j, k)は次式で計算できる。
Δφ1/2(i, j, k)＝arctan[ΔSi1/2(i, j, k)／ΔSr1/2(i, j, k)] ・・・・・（１６）
【００５３】
上記の本発明における３次元位相分布の計算に基づき、第２ステップ（102）から第４ステップ（104）の内容について説明する。
先ず、第２ステップ（102）においては、第１ステップ（101）で得られた第１のエコー信号(203,706)から計算した複素MR画像と第２のエコー信号(204,708)から計算した複素MR画像とから第１の複素演算MR画像を作成する。具体的には、第１ステップ（101）で得られた複素MR画像は式(3)で表されるもの、すなわち、S(i, j, k, 1)が第１のエコー信号から計算した複素MR画像、S(i, j, k, 2)が第２のエコー信号から計算した複素MR画像であり、第１の複素演算MR画像は式（５）で表されるものである。
【００５４】
次に、第３ステップ（103）では、各時相の第１の複素演算MR画像と、基準温度（加熱前温度）を反映した参照複素MR画像とから第２の複素演算MR画像を時系列に作成する。具体的には、加熱前の複素演算MR画像である式（７）のS1(i, j, k)と加熱後の複素演算画像である式（８）のS2(i, j, k)とから、式（９）または（１０）で表される第２の複素演算MR画像であるΔS1/2(i, j, k)を求めるものである。この具体的な例では、加熱前の第１の複素演算MR画像を基準にして、加熱後の第１の複素演算MR画像を時系列的に一定時間間隔（例えば、32秒間隔）で複数枚作成して、それらの各々の第１の複素演算MR画像と加熱前の第１の複素演算MR画像とから各々の第２の複素演算MR画像を作成することにより、加熱後の温度変化に対応する複数枚の時系列的な第２の複素演算MR画像が得られる。ここで、第２の複素演算MR画像を表す式（１０）からは、その位相成分を表す式（１１）が求まり、この式（１１）は加熱前後の温度変化に対応している（式（１２）、（１３）、（１a）参照）。
【００５５】
次に、第４ステップ（104）では、第２の複素演算MR画像のアークタンジェントを計算して、時系列位相分布を作成する。具体的には、式（１４）、（１５）、（１６）に基づいて、第３ステップ（103）で求めた時系列的な第２の複素演算MR画像から、各時相の位相成分Δφ1/2(i, j, k)を計算するものである。
【００５６】
ここで、逆三角関数arctanは、どのような値に対しても−πから＋πまでの値を与えるため、得られた位相差Δφ1-2(i, j, k)の画像はエリアシングを発生する。図３は、これを示したもので、図では見やすくするために２次元表示してある。以上の理由により、Δφ1-2(i, j, k)は、式（１）のθ(i, j, k)には直接対応していない。これを解決するために、次に述べる位相アンラップ処理を行う。
【００５７】
図１のフローチャートにおいて、次の第５ステップ（105）では、第４ステップ（104）で得られた３次元位相分布の計算結果について、３次元位相アンラップ処理を行う。ここで提案する位相アンラップ処理は、レーザ等によって被検体が局所的に加熱された場合、被検体内部の温度分布は連続的であるとの仮定に基づく。この仮定のもとに、３次元画像の被検体部分の非連続点を検出し、非連続点についてはエアリシングが発生していると判断して、位相値のアンラップ処理を行うものである。
【００５８】
以下、位相値のアンラップ処理について図３を参照しながら説明する。図３は本発明に用いる位相アンラップ処理を説明するための図である。位相アンラップ処理は、（１）基点の決定、（２）処理の走査方向の決定、（３）アンラップ処理の実行の手順で行われる。
【００５９】
(1)基点(i0, j0, k0)を決める。
図３の上図は被検体の位相図301を示したものである。この位相図301は２次元のもの（Δφ1-2(i, j, ko)）で、画素i、jについて表示されている。位相図301上には、加熱部分302、基点303、エリアシング304、305が表示されている。先ず、基点(i0, j0, k0)の決定方法としては、被検体内に存在し、かつ温度変化をしていない部位の１画素を選択する。温度変化をしている部位をおおまかに知るためには、温度上昇している領域の信号が低下することを利用する。すなわち、絶対値画像（式（７）、（８）参照）の加熱前後の信号差、
ΔSa1-2(i, j, k, 1)＝| Sa1(i, j, k, 1)−Sa2(i, j, k, 1)|・・・・・（１７）が、最大となる画素を温度変化があった部位と考え、基点から除外する。
【００６０】
基点303の決定については、元になる画像信号値Sa1(i, j, k, 1)およびSa2(i, j, k, 1)が大きいほうが、位相値の精度も高く、基点303として適している。すなわち、被検体内で温度変化（差分画像の絶対値の変化）が少なく、かつ元画像の信号量が大きい画素を基点303とするのが望ましい。
【００６１】
温度変化をしている部位をおおまかに知るための他の方法について図４により説明する。図４も図３の上図と同様に被検体の位相図401を示したものである。位相図401内には、温度変化した領域402、基点403が含まれている。この方法では、図４に示す如く撮影領域を数個（図４では９個）のブロックに分割し、各々のブロックで、被検体内のΔSa1-2(i, j, k, 1)の平均値Smeann（ｎ＝1〜ブロック数）を求め、Smeannが最大となるブロックn内を基点403の選定から除外する。この方法によれば、輝点ノイズなどに起因する不可避的な誤検出を避けることができる。
【００６２】
(2)処理の走査方法を決定する。
アンラップ処理の走査方向は最終的には全方向について行う。しかし、処理の順番は任意でよい。ここでは、先ず基点(i0, j0, k0)303を基準にしてiの正方向を走査方向とする。
【００６３】
(3)アンラップ処理を実行する。
基点(i0, j0, k0)303から出発し、走査方向に沿って順次、隣接画素間の位相差の差分を求め、位相値の差分が位相閾値thを越えていれば、アンラップ処理を行う。位相閾値thとしては、例えばπとする。具体的には、式（１８）に示すような演算を行う。
式(18)の演算内容を図３の下図を用いて説明する。先ず、図３の上図の被検体の位相図301のj＝j0線306上の位相差分Δφ1-2(i, j, k)が図３の下図の鋸歯状の実線Aである。この実線Aは、j＝j0線306上の点Cと点Dに対応する位置で位相値のジャンプをしている。
【００６５】
以下、順にiの正方向にi0+1、i0+2、i0+3、…について、アンラップ処理を繰り返す。被検体の点Eまでこの処理を行ったら、次に基点(i0, j0, k0)303を始点として、iの負方向にi0−1、i0−2、i0−3、…について、θ(i0, j0, k0)を初期値として同様の処理を行う。この結果、画素(i0, j0, k0)の線306についての、Δφ1-2(i, j0, k0)からθ(i, j0, k0)への変換が完了する。
【００６６】
ここで、位相アンラップ処理を正確に行うための前提条件として、「被検体の形状は、凸面で囲まれており、凹部はなく、従って、上記処理を行うにあたり、被検体領域がiの正または負方向に飛び飛びになることはない」がある。これは、通常の臨床利用ではこの条件が満たされていると考えられることから、実用上大きな制約となることはない。
【００６７】
次に、画素座標(i, j0, k0)を初期値として、処理の走査方向を、jもしくはkの方向とし、同様のアンラップ処理を行う。この結果、被検体内部について、３次元的な位相アンラップ処理が完了する。
【００６８】
位相アンラップ処理を更に安定にするためには、直交座標の３方向だけでなく、さらに斜方向、例えばi＝jと平行な方向にも同様なアンラップ処理を行うことが望ましい。
【００６９】
被検体外部の位相値は、背景ノイズから計算される位相値であり、ランダムな値である。これらの位相値は温度分布情報を含まないので、被検体外部の位相値θ(i, j, k)は0とする。
以上の位相アンラップ処理を行った結果、位相画像に存在したエリアシングが除去される。
【００７０】
図１のフローチャートにおいて、次の第６ステップ（106）では、第５ステップ（105）でアンラップ処理した３次元位相分布に基づいて３次元温度分布の計算を行う。すなわち、上記で計算したθ(i, j, k)に基づいて、被検体の温度変化の分布を、式（２）を３次元的表記した式（１２）を使って３次元的に求める。
ΔT(i, j, k)=θ(i, j, k)／(γ・B0・C・τ) ・・・・・（１９）
ここで、γ・B0は共鳴角周波数を与える。B0＝0.3Tの場合には、25.6×106・π[rad/s]である。τは撮影シーケンスのTE値（＝35ms）と等しい。
【００７１】
以上の如く、図１のフローチャートに示した手順に従うことにより、被検体の時系列温度分布を計算することができる。次に、この時系列温度分布の計算結果を表示する方法について説明する。
【００７２】
以下の実施例では、図１のフローチャートの第６ステップで求めた３次元温度分布の計算結果に基づき、表示のための３次元温度画像を計算し、これを表示画面に表示する。本発明においては、この表示に３次元投影法を採用する。温度分布の３次元投影法としては、ボリュームレンダリング法や最大投影値法などがあるが、ボリュームレンダリング法を用いるのが良い。何となれば、レーザによる被検体の加熱はレーザ照射部の中心部分の温度が最も高く、周囲に行くに従って温度が低くなるためである。
【００７３】
これに対し、最大投影値法を用いた場合、基本的には、加熱部分を３次元的に可視化できるものの、レーザファイバ先端部付近の閾値よりも高い高温度部が視覚的にわかりにくいと予想される。このため、本実施例では、温度画像をボリュームレンダリング法で可視化する。
【００７４】
ボリュームレンダリング法では、３次元空間に対してレイキャスティングを行う。この方法としては、例えば、front to back法、back to front法、並行投影法などが当業者に公知である。１画素に注目した場合、αの不透明度をもつ画素にCinの光が入って来たとき、画素の色をCとすると、画素を出て行く光Coutは式（２０）で表される。このCoutが撮影時のボクセル値に相当する。
Cout＝Cin(1−α)＋C・α ・・・・・（２０）
【００７５】
αは0で完全な透明、１で完全な不透明である。αは温度ΔT(i, j, k)に対応する。ΔT(i, j, k)＝0に対しては、α＝0とする。ΔT(i, j, k)の最大値は、例えばレーザ照射では300℃程度であるので、ΔT(i, j, k)＝300[℃]に対して、α＝1とする。すなわち、α＝ΔT(i, j, k)/300とする。このように設定することにより、レーザ照射点が不透明になり、その周辺に、楕円球状に広がる高温度領域を３次元的に画像化することができる。
【００７６】
温度計測のMR撮影を２回以上繰り返す場合には、n回目の撮影と１回目の撮影との間で、上記と同様の演算を行う。その結果、n回目の撮影と１回目の撮影の間に生じた温度変化の画像を３次元的に得ることができる。
【００７７】
さて、ボリュームレンダリング法は有用であるが、演算時間が長いため実用上使い勝手が若干悪いことは周知である。温度画像を連続的に、時系列的に、３次元表示したい場合に演算時間が長いことが問題となる。本実施例の撮影時間は32sまたは35.8sであるが、この数値は、例えばGrE型のマルチショットEPIシーケンスを導入することにより、1/4程度の時間、すなわち10s弱にまで短縮することが可能である。従って、３次元温度分布表示についても10s以内に完了することが重要な課題となる。
【００７８】
以下では、上記課題を解決するためのダイナミック３次元温度計測に適したボリュームレンダリングの改良手法について図５を参照しながら説明する。図５は本発明に用いるボリュームレンダリング処理の一実施例の一部を示す図である。ボリュームレンダリング法の演算時間を減らすには、元になる３次元温度画像の画素数を減らすことが本質的である。局所温度治療では、図５に示す如く、注目する温度変化部分501は、ほぼ球状に局在している。これは、レーザ照射のためのレーザファイバ502やRFトランスジューサが、侵襲的に患部は挿入され、その部位のみを加熱するためである。同様に、焦点超音波による加熱においても、局所のみの加熱となっている。
【００７９】
このように注目点が局在しているということは、血管像を表示するMRAの場合とは大きく異なる。温度画像のこの特徴を利用することによって、元の温度画像からボリュームレンダリングの対象となる画素のみを効率良く自動抽出することができる。
【００８０】
具体的には、温度画像ΔT(i, j, k)が表示される温度表示領域505のうちの、温度変化の最大点もしくは温度変化領域の重心点Q(imax, jmax, kmax)503を求め、この点Q503を中心としてユーザが設定する半径r0の球状領域504もしくは一辺の長さa0の立方体領域（図示せず）のみをボリュームレンダリングの対象領域とする。
【００８１】
球状領域504の半径r0は典型的には20mm程度、立方体領域の一辺の長さa0は40mm程度とする。撮影視野（FOV）が250mm×250mmの画像では、画素サイズは1.95mm(=250mm/128)なので、立方体領域の場合の一辺は20画素(＝40mm/1.95mm)である。また、ボリュームレンダリングに使われる画素数は、20×20×20となる。
【００８２】
ボリュームレンダリング対象領域504の一辺の大きさ（または半径の大きさ）は、撮影毎に変化させても良い。これは、レーザ加熱の場合、最初はごく小さい領域が加熱され、徐々に加熱領域が広がっていくことを考慮して、温度画像ΔT(i, j, k)の表示領域505のうちの、温度変化の最大点Q(imax, jmax, kmax)503を基準としてΔT(i, j, k)の空間的な変化をr(r＝|(i−imax, j−jmax, k−kmax)|)の関数として求める。ΔT(i, j, k)の値が最大変化点の値ΔT(imax, jmax, kmax)の半分になるような半径rを温度変化部位の広がりとして捕らえることができる。これをrdとする。
【００８３】
このときのボリュームレンダリング対象領域504の半径r0の値は、撮影毎に（すなわちボリュームレンダリングの元画像データの組毎に）、例えばr0＝3× rdとなるように変化させる。このように半径r0を設定することにより、ボリュームレンダリング対象の元画像は、温度変化領域の広がりとともに拡大していき、ボリュームレンダリング処理を行うのに必要十分な元画像データから構成され、ボリュームレンダリング処理の演算時間も最短となる。
【００８４】
上記の実施例では、レーザによる加熱時の被検体の温度変化を３次元的に可視化する例を示したが、他の実施例として、例えば、RFアブレーション時の被検体の温度変化を画像化することもできる。また、例えばクライオサージェリーでは、−40℃まで患部を冷却した場合、冷凍領域の直径は20mm程度に達するので、このような低温治療時の温度モニタにも本発明は利用することができる。この場合、ボリュームレンダリングにおいて、温度が下がるほどαの値を大きくするように変換式の符号を反転する必要がある。
【００８５】
【発明の効果】
以上説明した如く、本発明によれば、MRI装置を用いて被検体内の温度分布を計測する場合に、被検体の動きがあっても、被検体内の温度変化を安定に計算することができる。
【００８６】
また、従来生体の加熱時の温度分布を計測するためには、加熱前と加熱後に撮影し、２画像の位相差分を求める必要があり、この間に被検体が動くと計算結果に誤差を発生しやすかった。しかし、本発明によれば、加熱前後の時間が5秒ないし10秒必要となる場合や加熱に伴う痛みがある場合などに、被検体が動くことがあったり、呼吸や心拍、腸の蠕動などに起因する動きが発生しても、上記の画像差分時に臓器の境界部分に生じていたアーチファクトを低減することができる。
【００８７】
また、本発明によれば、従来報告されているオフライン処理の体動抑制法と異なり、オンラインで体動アーチファクトを低減することができる。その結果、温度計測のように、準リアルタイムで結果を見ながら手術を行うIV-MRIに好適な画像が得られる。また、伸縮性の幅広の帯を使う方法と異なり、手術部位を覆ってしまうこともないのでIV-MRIに適している。
【００８８】
また、従来はレーザ照射などで得られる加熱部分を３次元で計測し、３次元で表示する手法は検討されていなかったが、本発明によれば、レーザ照射時の局所的に急激な位相変化がある場合に、立体構造を安定かつ正確に表示することができる。また、位相から血管像を求めるPC（Phase contrast）法のように撮影シーケンスのパラメータを、最適にセットする必要もなく、３次元温度計測が可能である。
【００８９】
また、本発明によれば、加熱部分の時系列的な変化を準リアルタイムで３次元的に可視化できるMRI装置を提供でき、加熱治療時の３次元的な温度変化の様子を直ちにチェックすることができる。
【図面の簡単な説明】
【図１】本発明の時系列温度計測方法の一実施例のフローチャート。
【図２】本発明に用いる撮影シーケンスの一実施例。
【図３】本発明に用いるエリアシング処理を説明するための図。
【図４】本発明に用いる位相アンラップ処理を説明するための他の図。
【図５】本発明に用いるボリュームレンダリング処理の一実施例の一部を示す図。
【図６】本発明を実施するのに用いるMRI装置の構成例を示すブロック図。
【図７】本発明に用いる撮影シーケンスの他の実施例。
【図８】ペネトレーション効果を示す模式図。
【符号の説明】
101…第１ステップ
102…第２ステップ
103…第３ステップ
104…第４ステップ
105…第５ステップ
106…第６ステップ
201, 701…90度RFパルス
202…180度RFパルス
203, 706…第１のエコーパルス
204, 708…第２のエコーパルス
205…TE／2
206…τ／2
207, 208, 705, 707…読み出し傾斜磁場
209, 702…スライス傾斜磁場
210, 704…位相エンコード傾斜磁場
301, 401…位相図
302, 402…加熱部分
303, 403…基点
304, 305…エリアシング
306…j＝j0線
501…温度変化部分
502…レーザファイバ
503…温度変化の最大点もしくは温度変化領域の重心点Q
504…球状領域(ボリユームレンダリング対象領域)
505…温度表示領域
600…計測空間
601…被検体
602…磁石
603…傾斜磁場コイル
604…RF照射コイル
605…RF受信コイル
606…信号検出部
607…信号処理部
608…表示部
609…傾斜磁場電源
610…RF送信部
611…制御部
612…ベッド
703…スライスエンコード傾斜磁場
709…第１のエコー時間
710…第２のエコー時間
801, 804…RF磁場
802…被検体
803…渦電流
Claims (4)
- 被検体から計測されたエコー時間の異なる第１のエコー信号と第２のエコー信号とを用いて、前記被検体の温度変化を計測する方法であって、
（１）前記第１のエコー信号と前記第２のエコー信号とを用いて、第１の複素画像を作成するステップと、
（２）前記第１の複素画像と基準温度を反映した参照複素画像とから、第２の複素画像を作成するステップと、
（３）前記第２の複素画像から位相分布を作成するステップと、
（４）前記位相分布に基づき、前記被検体の温度変化分布を取得するステップと、
を含む温度計測方法であって、
前記第１のエコー信号はスピンエコー時間より前に計測され、前記第２のエコー信号は前記スピンエコー時間の後に計測されたエコー信号であることを特徴とする温度計測方法。 - 請求項１記載の温度計測方法において、
前記（１）〜（４）のステップを繰り返して、前記被検体の温度変化分布を時系列に取得することを特徴とする温度計測方法。 - 請求項１又は２に記載の温度計測方法において、
前記基準温度は、前記被検体を加熱する前の温度であり、前記第２の複素画像は、前記被検体を加熱した状態の画像であることを特徴とする温度計測方法。 - 被検体からエコー時間の異なる第１のエコー信号と第２のエコー信号の計測を制御する計測制御手段と、
前記第１のエコー信号と前記第２のエコー信号とを用いて、前記被検体の温度変化を取得する演算処理手段と、
を備え、
前記演算処理手段は、
前記第１のエコー信号と前記第２のエコー信号とを用いて、第１の複素画像を作成し、
前記第１の複素画像と基準温度を反映した参照複素画像とから、第２の複素画像を作成し、
前記第２の複素画像から位相分布を作成し、
前記位相分布に基づき、被検体の温度変化分布を取得する磁気共鳴イメージング装置であって、
前記第１のエコー信号はスピンエコー時間より前に計測され、前記第２のエコー信号は前記スピンエコー時間の後に計測されたエコー信号であることを特徴とする磁気共鳴イメージング装置。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP11549999A JP4250255B2 (ja) | 1999-04-22 | 1999-04-22 | 温度計測方法及び磁気共鳴イメージング装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP11549999A JP4250255B2 (ja) | 1999-04-22 | 1999-04-22 | 温度計測方法及び磁気共鳴イメージング装置 |
Publications (3)
| Publication Number | Publication Date |
|---|---|
| JP2000300535A JP2000300535A (ja) | 2000-10-31 |
| JP2000300535A5 JP2000300535A5 (ja) | 2006-08-10 |
| JP4250255B2 true JP4250255B2 (ja) | 2009-04-08 |
Family
ID=14664033
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP11549999A Expired - Fee Related JP4250255B2 (ja) | 1999-04-22 | 1999-04-22 | 温度計測方法及び磁気共鳴イメージング装置 |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP4250255B2 (ja) |
Families Citing this family (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1993628B (zh) * | 2004-08-02 | 2011-05-04 | 皇家飞利浦电子股份有限公司 | 一种磁共振成像系统和计算温度差值图的方法 |
| CN102772207B (zh) * | 2011-05-12 | 2015-05-13 | 上海联影医疗科技有限公司 | 磁共振成像方法和装置 |
| CN104981709B (zh) * | 2012-05-02 | 2018-02-02 | 皇家飞利浦有限公司 | 成像热测量 |
| DE102012221463B4 (de) | 2012-11-23 | 2023-09-28 | Siemens Healthcare Gmbh | Verfahren und Vorrichtung zur schnellen Erkennung von Bereichen geänderter Temperatur |
| JP6595393B2 (ja) * | 2016-04-04 | 2019-10-23 | 株式会社日立製作所 | 磁気共鳴イメージング装置、及び、画像処理方法 |
| EP3378426A1 (en) * | 2017-03-20 | 2018-09-26 | Koninklijke Philips N.V. | Locating ablated tissues using electric properties tomography |
-
1999
- 1999-04-22 JP JP11549999A patent/JP4250255B2/ja not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| JP2000300535A (ja) | 2000-10-31 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US7542793B2 (en) | MR-guided breast tumor ablation and temperature imaging system | |
| US7725157B2 (en) | System and method for interventional procedures using MRI | |
| EP2499508B1 (en) | Correcting artifacts caused by motion of a medical device and/or patient motion in magnetic resonance thermometry | |
| KR101461099B1 (ko) | 자기공명영상장치 및 기능적 영상획득방법 | |
| RU2580189C2 (ru) | Интервенционная мр-томография с компенсацией движения | |
| US20050065429A1 (en) | Method for three plane interleaved acquisition for three dimensional temperature monitoring with MRI | |
| JPWO2002022012A1 (ja) | 磁気共鳴イメージング装置 | |
| JP3976845B2 (ja) | 磁気共鳴イメージング装置 | |
| JP4250255B2 (ja) | 温度計測方法及び磁気共鳴イメージング装置 | |
| JP4250251B2 (ja) | 磁気共鳴画像診断装置 | |
| JP4349647B2 (ja) | 磁気共鳴イメージング装置 | |
| JP2000300536A (ja) | Ｍｒｉ装置を用いた３次元温度計測方法 | |
| JP3817383B2 (ja) | Ｍｒｉ装置 | |
| KR20110044161A (ko) | 키홀 방법을 이용한 자기공명 온도 감시 영상 획득 방법 | |
| JP4177165B2 (ja) | Ｍｒｉ装置 | |
| JP4219032B2 (ja) | 磁気共鳴画像診断装置 | |
| JP5371620B2 (ja) | 核磁気共鳴イメージング装置 | |
| JP4371510B2 (ja) | 磁気共鳴イメージング装置 | |
| JP4118119B2 (ja) | 磁気共鳴イメージング装置 | |
| JP7357516B2 (ja) | 磁気共鳴撮像装置、及び、その制御方法 | |
| KR101480413B1 (ko) | B1 정보 획득 방법 및 장치 | |
| JP6912341B2 (ja) | 磁気共鳴イメージング装置、それを用いたデバイス位置検出方法、および、画像ガイド下インターベンション支援装置 | |
| US20240353514A1 (en) | Systems and methods for magnetic resonance based skull thermometry | |
| JP4074513B2 (ja) | 磁気共鳴イメージング装置 | |
| Matsuoka et al. | Integrated MR-laparoscopy system with respiratory synchronization for minimally invasive liver surgery |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20060419 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20060425 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20060517 |
|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20060419 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20060517 |
|
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20080213 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20080225 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20080410 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20080821 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20081014 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20081110 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20081210 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20090113 |
|
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20090119 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120123 Year of fee payment: 3 |
|
| R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130123 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20140123 Year of fee payment: 5 |
|
| LAPS | Cancellation because of no payment of annual fees |