JP2004081657A

JP2004081657A - Method for extracting fibrous image, image processing device, and magnetic resonance imaging systems

Info

Publication number: JP2004081657A
Application number: JP2002248285A
Authority: JP
Inventors: Kiyoko Abe; 阿部　清子; Hiroyuki Kabasawa; 椛沢　宏之
Original assignee: GE Medical Systems Global Technology Co LLC
Current assignee: GE Medical Systems Global Technology Co LLC
Priority date: 2002-08-28
Filing date: 2002-08-28
Publication date: 2004-03-18

Abstract

<P>PROBLEM TO BE SOLVED: To realize a method for extracting fibrous images which enables an automatic and rapid extraction of nerve fiber images, an image processing device, and a magnetic resonance imaging system. <P>SOLUTION: A diffusion tensor matrix is calculated from information on diffusion tensor images and T<SB>2</SB>weighted images (Step S603). Eigenvalues and eigenvectors of the diffusion tensor matrix are obtained by diagonalyzing the matrix (Step S604). A fractional anisotropy (FA) value which is an index of anisotropy of the diffusion coefficient is calculated from the eigenvalues (Step S605). A region of the nerve fibers is extracted from the FA value (Step S606). Since the pixels of the nerve fibers in the direction of the maximum eigenvector showing the direction of the nerve fibers are extracted sequentially from the region of the nerve fibers by using a region extension means 250 (Step S607), the extraction of the information on the nerve fiber images is realized in a short time without individual differences. <P>COPYRIGHT: (C)2004,JPO

Description

【０００１】
【発明の属する技術分野】
この発明は、ディフュージョンテンソル（ｄｉｆｆｕｓｉｏｎ　ｔｅｎｓｏｒ）画像情報を用いて、撮像部位から抽出される繊維状画像抽出方法、画像処理装置および磁気共鳴撮像システムに関する。
【０００２】
【従来の技術】
近年、磁気共鳴撮像装置の分野においては、従来の被検体の形態的な画像情報に加えて、被検体の機能を画像化したファンクショナル（ｆｕｎｃｔｉｏｎａｌ）画像が臨床的に用いられつつある。このファンクショナル画像の中で、拡散強調画像は、組織内の拡散係数をディフュージョンテンソル画像から算出し、画像化する方法である。また、この拡散係数の異方性は、脳あるいは脊髄の白質神経繊維においては、繊維の走行方向を示していることが知られている（阿部修等：ＭＲＩ：拡散画像；画像診断、Ｖｏｌ１９、Ｎｏ．１０、１９９９参照）。
【０００３】
そこで、この拡散係数の異方性の大小から、神経繊維が存在する領域を選択し、さらにこの選択領域から、神経繊維の走行方向を示す神経繊維画像を抽出することが行われている。
【０００４】
【発明が解決しようとする課題】
しかしながら、上記従来技術によれば、神経繊維画像のみの抽出には、多大な時間と専門的な知識とが必要である。すなわち、神経繊維画像を選択領域から抽出する作業は、医者が手作業により、解剖学的知識に基づいて、画素ごとに取捨選択を行う必要がある。
【０００５】
特に、医者が手作業で、神経繊維画像のみの抽出を行う際には、解剖学的知識に依存する個人差が生じると共に、多量の画像情報を取捨選択するのに多大な時間を要し、この長時間の作業に起因する人為的な間違いが生じる要因となっている。
【０００６】
これらのことから、神経繊維画像のみの抽出を、自動で高速に行える繊維状画像抽出方法、画像処理装置および磁気共鳴撮像システムをいかに実現するかが重要となる。
【０００７】
この発明は、上述した従来技術による課題を解決するためになされたものであり、神経繊維画像のみの抽出を、自動で高速に行える繊維状画像抽出方法、画像処理装置および磁気共鳴撮像システムを提供することを目的とする。
【０００８】
【課題を解決するための手段】
上述した課題を解決し、目的を達成するために、第１の観点の発明にかかる繊維状画像抽出方法は、イメージングパルスシーケンスを実行して被検体の同一撮像部位のディフュージョンテンソル画像情報およびＴ_２強調画像情報を取得し、前記ディフュージョンテンソル画像情報および前記Ｔ_２強調画像情報から画素ごとの拡散テンソル行列を算出し、前記拡散テンソル行列の固有値および固有ベクトルを算出し、前記固有値から前記画素位置における拡散係数の異方性の指標値を算出し、前記指標値に基づいて、前記撮像部位から異方性領域を抽出し、前記異方性領域内の前記指標値が閾値を越える領域に始点画素を選定し、前記画素ごとの前記固有ベクトルの中から、最大の前記固有値を有する最大固有ベクトルを選択し、前記始点画素から、前記最大固有ベクトルのベクトル情報を用い、繊維状画像を順次抽出することを特徴とする。
【０００９】
この第１の観点による発明によれば、イメージングパルスシーケンスを実行して被検体の同一撮像部位のディフュージョンテンソル画像情報およびＴ_２強調画像情報を取得し、このディフュージョンテンソル画像情報およびＴ_２強調画像情報から画素ごとの拡散テンソル行列を算出し、この拡散テンソル行列の固有値および固有ベクトルを算出し、この固有値から画素位置における拡散係数の異方性の指標値を算出し、この指標値に基づいて、撮像部位から異方性領域を抽出し、この異方性領域内の指標値が閾値を越える領域に始点画素を選定し、画素ごとの固有ベクトルの中から、最大の固有値を有する最大固有ベクトルを選択し、始点画素から、最大固有ベクトルのベクトル情報を用い、繊維状画像を順次抽出することとしているので、神経繊維の進行方向を示す最大固有ベクトルを算出し、このベクトルの方向に沿った神経繊維画素を順次，自動で抽出することができる。
【００１０】
また、第２の観点の発明にかかる画像処理装置は、被検体の同一撮像部位のディフュージョンテンソル画像情報およびＴ_２強調画像情報を記憶する記憶手段と、前記ディフュージョンテンソル画像情報および前記Ｔ_２強調画像情報から画素ごとの拡散テンソル行列を算出するテンソル行列算出手段と、前記拡散テンソル行列の固有値および固有ベクトルを算出する行列対角化手段と、前記固有値から前記画素位置における拡散係数の異方性の指標値を算出する指標算出手段と、前記指標値に基づいて、前記撮像部位から高い異方性領域を抽出する閾値手段と、を備える画像処理装置であって、前記異方性領域内の前記指標値が閾値を越える領域に始点画素を選定し、前記画素ごとの前記固有ベクトルの中から、最大の前記固有値を有する最大固有ベクトルを選択し、前記始点画素から、前記最大固有ベクトルのベクトル情報を用いて、繊維状画像を抽出する領域拡張手段と、を備えることを特徴とする。
【００１１】
この第２の観点の発明によれば、領域拡張手段を用いて、異方性領域内の指標値が閾値を越える領域に始点画素を選定し、画素ごとの固有ベクトルの中から、最大の固有値を有する最大固有ベクトルを選択し、始点画素から、最大固有ベクトルのベクトル情報を用いて、繊維状画像を抽出することとしているので、神経繊維の進行方向を示す最大固有ベクトルを算出し、その方向の画素を領域拡張手段により抽出して行くことにより、神経繊維画像を自動抽出することができる。
【００１２】
また、第３の観点の発明にかかる画像処理装置は、前記指標値として、前記固有値から算出されるフラクショナルアニソトロピー値を用いることを特徴とする。
【００１３】
この第３の観点の発明によれば、指標値として、固有値から算出されるフラクショナルアニソトロピー値を用いることとしているので、拡散係数の異方性を定量的に評価することができる。
【００１４】
また、第４の観点の発明にかかる画像処理装置は、前記閾値手段により、前記フラクショナルアニソトロピー値が閾値を越える画素領域を抽出することを特徴とする。
【００１５】
この第４の観点の発明によれば、閾値手段は、フラクショナルアニソトロピー値が閾値を越える画素領域を抽出することとしているので、異方性が高い神経繊維領域を、フラクショナルアニソトロピー値により抽出することができる。
【００１６】
また、第５の観点の発明にかかる画像処理装置は、前記テンソル行列算出手段が、最小二乗法を用いて最適化されることを特徴とする。
この第５の観点の発明によれば、テンソル行列算出手段は、最小二乗法を用いて最適化されることとしているので、確度の高い拡散テンソル行列を求めることができる。
【００１７】
また、第６の観点の発明にかかる画像処理装置は、前記領域拡張手段が、マスク領域の中心画素と周辺画素とのつながりを示す結合指数を、前記ベクトル情報から前記周辺画素ごとに算出し、前記結合指数が閾値を越えない前記周辺画素を新たな中心画素とする新たなマスク領域を設定し、順次領域拡張を行うことを特徴とする。
【００１８】
この第６の観点の発明によれば、領域拡張手段が、マスク領域の中心画素と周辺画素とのつながりを示す結合指数を、ベクトル情報から周辺画素ごとに算出し、この結合指数が閾値を越えない周辺画素を新たな中心画素とする新たなマスク領域を設定し、順次領域拡張を行うこととしているので、中心画素の神経繊維と同一の神経繊維からなる周辺画素を選択し、領域拡張を用いて、神経繊維全体を抽出することができる。
【００１９】
また、第７の観点の発明にかかる画像処理装置は、前記結合指数に、前記中心画素および前記周辺画素の前記ベクトル情報の内積から算出される方向余弦の角度を用いることを特徴とする。
【００２０】
この第７の観点の発明によれば、結合指数に、中心画素および周辺画素のベクトル情報の内積から算出される方向余弦の角度を用いることとしているので、最大固有ベクトルの向きが概ね同一の画素を、結合指数から選択することができる。
【００２１】
また、第８の観点の発明にかかる磁気共鳴撮像システムは、静磁場を形成し、勾配磁場を形成し、ＲＦ信号を送信あるいは受信し、前記勾配磁場および前記ＲＦ信号を制御して、被検体の同一撮像部位のディフュージョンテンソル画像情報およびＴ_２強調画像情報を取得する磁気共鳴撮像装置と、前記磁気共鳴撮像装置の前記ディフュージョンテンソル画像情報およびＴ_２強調画像情報を転送する転送手段と、前記転送手段により転送されたディフュージョンテンソル画像情報およびＴ_２強調画像情報に画像処理を行う画像処理装置と、を備える磁気共鳴撮像システムであって、前記画像処理装置は、前記ディフュージョンテンソル画像情報および前記Ｔ_２強調画像情報から画素ごとの拡散テンソル行列を算出するテンソル行列算出手段と、前記拡散テンソル行列の固有値および固有ベクトルを算出する行列対角化手段と、前記固有値から前記画素位置における拡散係数の異方性の指標値を算出する指標算出手段と、前記指標値に基づいて、前記撮像部位から異方性領域を抽出する閾値手段と、前記異方性領域内の前記指標値が閾値を越える領域に始点画素を選定し、前記画素ごとの前記固有ベクトルの中から、最大の固有値を有する最大固有ベクトルを選択し、前記始点画素から、前記最大固有ベクトルのベクトル情報を用いて、繊維状画像を抽出する領域拡張手段と、を備えることを特徴とする。
【００２２】
この第８の観点の発明によれば、画像処理装置は、テンソル行列算出手段により、ディフュージョンテンソル画像情報およびＴ_２強調画像情報から画素ごとの拡散テンソル行列を算出し、行列対角化手段により、この拡散テンソル行列の固有値および固有ベクトルを算出し、指標算出手段により、この固有値から画素位置における拡散係数の異方性の指標値を算出し、この指標値に基づいて、閾値手段により、撮像部位から異方性領域を抽出し、領域拡張手段により、この異方性領域内の指標値が閾値を越える領域に始点画素を選定し、画素ごとの固有ベクトルの中から、最大の固有値を有する最大固有ベクトルを選択し、始点画素から、最大固有ベクトルのベクトル情報を用いて、繊維状画像を抽出することとしているので、神経繊維の進行方向を示す最大固有ベクトルを算出し、その方向の画素を領域拡張手段により抽出して行くことにより、神経繊維画像を自動抽出することができる。
【００２３】
【発明の実施の形態】
以下に添付図面を参照して、この発明にかかる繊維状画像抽出方法、画像処理装置および磁気共鳴撮像システムの好適な実施の形態について説明する。なお、これにより本発明が限定されるものではない。
【００２４】
まず、本実施の形態にかかる磁気共鳴撮像システムの全体構成について説明する。図１に磁気共鳴撮像システムのブロック（ｂｌｏｃｋ）図を示す。同図に示すように、本システムはマグネット（ｍａｇｎｅｔ）部１００、キャビネット（ｃａｂｉｎｅｔ）部２０、操作コンソール（ｃｏｎｓｏｌｅ）部３０、画像処理装置４０からなる。
【００２５】
マグネット部１００は、主磁場コイル（ｃｏｉｌ）部１０２、勾配コイル部１０６およびＲＦ（ｒａｄｉｏ　ｆｒｅｑｕｅｎｃｙ）コイル部１０８を有する。これら各コイル部は概ね円筒形の形状を有し、互いに同軸的に配置されている。また、マグネット部１００の概ね円柱状の撮像空間（ボア：ｂｏｒｅ）の中心部には、クレードル（ｃｒａｄｌｅ）１２０上に横臥状態で撮像される被検体１が配置される。
【００２６】
主磁場コイル部１０２は、マグネット部１００の内部空間に静磁場を形成する。静磁場の方向は概ね同軸的に配置された円筒状コイルの軸方向に平行である。すなわち、いわゆる水平磁場を形成する。主磁場コイル部１０２は、例えば超伝導コイルを用いて構成される。なお、超伝導コイルに限らず常伝導コイル等を用いて構成してもよい。
【００２７】
勾配コイル部１０６は、図示しない互いに独立な３系統の勾配コイル、すなわち、ｘ、ｙ、ｚ座標軸方向の静磁場強度に線形勾配を持たせるための３つの勾配コイルからなる。ここで、スライス（ｓｌｉｃｅ）軸、位相エンコード（ｅｎｃｏｄｅ）軸および周波数エンコード軸は、撮像断面を決定する、互いに直交する座標軸で相互間の垂直性を保ったままｘ、ｙ、ｚ座標軸に関して任意の傾きを持たせることも可能である。
【００２８】
スライス軸方向の勾配磁場をスライス勾配磁場、位相エンコード軸方向の勾配磁場を位相エンコード勾配磁場、周波数エンコード軸方向の勾配磁場を周波数エンコード勾配磁場ともいう。以下、勾配磁場を単に勾配ともいう。
【００２９】
ＲＦコイル部１０８は、静磁場空間に、被検体１の体内のスピン（ｓｐｉｎ）を励起するための高周波磁場を形成する。以下、高周波磁場を形成することをＲＦ励起信号の送信ともいう。また、ＲＦ励起信号をＲＦパルス（ｐｕｌｓｅ）ともいう。ＲＦコイル部１０８は、ＲＦ励起信号の送信と共に、励起されたスピンが生じる電磁波すなわち磁気共鳴信号を受信する。
【００３０】
ＲＦコイル部１０８は、送信用のコイルおよび受信用のコイルを有する。送信用のコイルおよび受信用のコイルは、同じコイルを兼用するかあるいはそれぞれ専用のコイルを用いる。
【００３１】
キャビネット部２０は、スキャンコントローラ（ｓｃａｎ　ｃｏｎｔｒｏｌｌｅｒ）部１６０、勾配駆動部１３０、ＲＦ駆動部１４０、データ（ｄａｔａ）収集部１５０、データ処理部１５５、記憶部１６６からなる。
【００３２】
スキャンコントローラ部１６０は、勾配駆動部１３０、ＲＦ駆動部１４０、データ収集部１５０、データ処理部１５５を、制御する。
勾配駆動部１３０には、勾配コイル部１０６が接続されている。勾配駆動部１３０は、勾配コイル部１０６に駆動信号を与えて勾配磁場を発生させる。勾配駆動部１３０は、勾配コイル部１０６における３系統の勾配コイルに対応して、図示しない３系統の駆動回路を有する。
【００３３】
ＲＦ駆動部１４０には、ＲＦコイル部１０８が接続されている。ＲＦ駆動部１４０は、ＲＦコイル部１０８に駆動信号を与えてＲＦパルスを送信し、被検体１の体内のスピンを励起する。
【００３４】
また、データ収集部１５０も、ＲＦコイル部１０８に接続されている。データ収集部１５０は、ＲＦコイル部１０８が受信した受信信号をサンプリング（ｓａｍｐｌｉｎｇ）によって取り込み、それをディジタルデータ（ｄｉｇｉｔａｌ　ｄａｔａ）として収集する。データ収集部１５０の出力側は、記憶部１６６に接続され、これらディジタルデータは、記憶部１６６に、生データとして、蓄積、保存される。
【００３５】
スキャンコントローラ部１６０は、例えばコンピュータ（ｃｏｍｐｕｔｅｒ）等を用いて構成される。スキャンコントローラ部１６０は、図示しないメモリ（ｍｅｍｏｒｙ）を有する。メモリはスキャンコントローラ部１６０用の制御プログラム（ｐｒｏｇｒａｍ）であるパルスシーケンス（ｐｕｌｓｅ　ｓｅｑｕｅｎｃｅ）および各種のデータを記憶している。スキャンコントローラ部１６０の機能は、コンピュータがメモリに記憶されたパルスシーケンスを実行することにより実現される。
【００３６】
記憶部１６６には、ディジタル化された受信信号が、生データとして蓄積されており、この受信信号は、撮像の終了とともに、データ処理部１５５に転送され画像再構成が行われる。データ処理部１５５は、コンピュータ等を用いて構成され、データ収集部１５０が収集した磁気共鳴信号である、Ｋ空間上の生データを２次元フーリエ変換により画像再構成を行い、画像データを生成する。そして、この画像データを記憶部１６６に転送し、保存する。
【００３７】
操作コンソール部３０の制御部１７０は、キャビネット部２０のスキャンコントローラ部１６０およびデータ処理部１５５の上位にあって、これらを制御する。また、これらにお動作は、制御部１７０が記憶部１７６に記憶されたプログラムを実行することによりを実現される。
【００３８】
操作コンソール部３０は、制御部１７０、記憶部１７６、表示部１８０および操作部１９０からなる。記憶部１７６は、キャビネット部２０の記憶部１６６と、図示しない通信手段により接続され、記憶部１６６の画像データを取得する。表示部１８０は、グラフィックディスプレー（ｇｒａｐｈｉｃ　ｄｉｓｐｌａｙ）等で構成される。操作部１９０はポインティングデバイス（ｐｏｉｎｔｉｎｇｄｅｖｉｃｅ）を備えたキーボード（ｋｅｙｂｏａｒｄ）等で構成される。
【００３９】
表示部１８０は、記憶部１７６の画像データおよび各種の情報を表示する。操作部１９０は、オペレータ（ｏｐｅｒａｔｏｒ）によって操作され、各種の指令や情報等を制御部１７０に入力する。オペレータは、表示部１８０および操作部１９０を通じてインタラクティブ（ｉｎｔｅｒａｃｔｉｖｅ）に本装置を操作する。
【００４０】
画像処理装置４０は、記憶部１７７、画像処理部１７１、表示部１８１および操作部１９１を含んでいる。記憶部１７７は、操作コンソール部３０の記憶部１７６と、図示しない通信手段により接続され、記憶部１７６の画像データを取得する。画像処理部１７１は、記憶部１７７の画像データに画像処理を施し、表示部１８１に表示する。表示部１８１は、グラフィックディスプレー等で構成される。操作部１９１はポインティングデバイスを備えたキーボード等で構成される。
【００４１】
表示部１８１は、記憶部１７７の画像データ、画像処理データおよび各種の情報を表示する。操作部１９１は、オペレータによって操作され、各種の指令や情報等を画像処理部１７１に入力する。オペレータは、表示部１８１および操作部１９１を通じてインタラクティブに本装置を操作する。
【００４２】
ここで、画像処理装置４０の機能ブロック図および動作を説明する前に、エコープラナー（ｅｃｈｏ　ｐｌａｎａｒ）法を用いたディフュージョンテンソル画像の取得について、図２および３を用いて説明する。まず、ディフュージョンテンソル画像は、被検体１中に存在する水プロトン（ｐｒｏｔｏｎ）の拡散係数の分布を算出する基となるデータの画像情報である。ここで、被検体１中の拡散係数は、同一画素位置でも方向により異なる値を有し、行列を用いたテンソル（Ｔｅｎｓｏｒ）量により記述される。
【００４３】
オペレータは、操作部１９０から、被検体１の撮像面および撮像のスライス枚数を設定する。図２（Ａ）は、被検体１の頭部をＮ枚のチルト（ｔｉｌｔ）されたアキシャル（ａｘｉａｌ）撮像面にてスライスする例をサジタル（ｓａｇｉｔｔａｌ）断面から示した図である。このスライス面およびスライス枚数は、撮像する神経繊維領域の範囲および神経繊維の向きを考慮して決定される。
【００４４】
その後、運動検出勾配磁場（以下、ＭＰＧ勾配：ｍｏｔｉｏｎ　ｐｒｏｂｉｎｇ　ｇｒａｄｉｅｎｔ　と称する）が付加されたスピンエコーを用いたエコープラナー法を用いて、スライスごとの画像情報を取得する。このエコープラナー法のパルスシーケンスの例を図２（Ｂ）に示した。図２（Ｂ）は、横軸は共通の時間軸を有し、互いに独立するスライス軸、位相エンコード軸および周波数エンコード軸の３軸方向の動作を示すタイムチャート（ｔｉｍｅ　ｃｈａｒｔ）である。
【００４５】
図２（Ｂ）のＲＦパルスは、９０度パルスにより、スライス勾配で指定された被検体１の目的とするスライス面の励起を行い、さらに１８０度パルスにより励起プロトンの位相反転を行い、スピンエコーを取り出す。
【００４６】
また、位相エンコード勾配は、位相エンコード軸方向の励起プロトンに位相変化を加え、この位相変化に対応する周波数成分情報を取り出す。また、周波数エンコード勾配は、読み出し勾配とも呼ばれ、周波数エンコード軸方向に勾配を加え、この勾配の印可中に観測されるスピンエコーの周波数成分から、被検体１内部の位置が特定される。なお、図２（Ａ）に示した様な、被検体１の体軸方向と角度をなすスライス面が選択される際には、位相エンコード勾配、周波数エンコード勾配およびスライス勾配は、ｘ、ｙ、ｚ座標軸方向と一致せず、相互間の垂直性を保ったまま、ベクトルの座標変換と同様の変換により、ｘ、ｙ、ｚ座標軸に分配される。
【００４７】
また、図２（Ｂ）のＭＰＧ勾配は、１８０度パルスを挟んで対称位置に配設された、同一の大きさおよび時間間隔を有する勾配である。ＭＰＧ勾配は、図２（Ｂ）では、スライス軸、位相エンコード軸および周波数エンコード軸とは別の軸に示されているが、実際のスキャンでは、スライス軸、位相エンコード軸および周波数エンコード軸の少なくとも１つに大きさを変化させて付加される。これにより、拡散係数に寄与する、規則性を持たないランダム（ｒａｎｄｏｍ）な動きをする水プロトンのエコー信号は、位相の拡散により減少する。さらに、ＭＰＧ勾配の大きさおよび付加される軸方向を変化させることにより拡散係数の異方性情報も取得される。
【００４８】
ここで、ｘ座標軸方向に印可されるＭＰＧ勾配の大きさをｂｘ、ｙ座標軸方向に印可されるＭＰＧ勾配の大きさをｂｙ、ｚ座標軸方向に印可されるＭＰＧ勾配の大きさをｂｚとする。この勾配は、スライス軸、位相エンコード軸、周波数エンコード軸、あるいはｘ、ｙ、ｚ座標軸の同一時間軸上に配設される。そして、各ｘ、ｙ、ｚ座標軸のＭＰＧ勾配の大きさの組み合わせ（ｂｘ、ｂｙ、ｂｚ）ごとに、Ｎ枚からなるスライス面の生データが取得される。
【００４９】
図３に、異なるＭＰＧ勾配のｍ個の組み合わせの撮像から取得されるディフュージョンテンソル画像情報２００の例を模式的に示した。ＭＰＧ勾配の組み合わせごとに、図２（Ｂ）のパルスシーケンスで取得された、スライス面数であるＮ枚の生データを、データ処理部１５５で画像再構成した画像データである。なお、後述する拡散テンソル行列を算出するために、ｍ＞６である必要がある。また、異なるＭＰＧ勾配の組み合わせは、添字ｍにより示した。
【００５０】
つづいて、画像処理装置４０の機能構成について説明する。図４は、画像処理装置４０の機能ブロック図を示した図である。画像処理装置４０は、画像情報であるディフュージョンテンソル画像情報２００、Ｔ_２強調画像情報２０２、演算手段であるテンソル行列算出手段２１０、行列対角化手段２２０、指標算出手段２３０、領域設定手段である閾値手段２４０、領域拡張手段２５０、および表示手段であるボリュームレンダリング（ｖｏｌｕｍｅ　ｒｅｎｄｅｒｉｎｇ）手段２６０等からなる。
【００５１】
ディフュージョンテンソル画像情報２００は、キャビネット部２０の記憶部１６６から、操作コンソール部３０の記憶部１７６に転送され、さらに画像処理装置４０の記憶部１７７に転送される。このディフュージョンテンソル画像情報２００は、図３に示した様に、ＭＰＧ勾配の大きさの組み合わせの数、例えばｍ個からなる組み合わせごとの、スライス枚数、例えばＮ枚の、画像データから構成される。
【００５２】
Ｔ_２強調画像情報２０２は、ディフュージョンテンソル画像情報２００のＭＰＧ勾配の大きさが、すべての座標軸で零となる場合に、図２（Ｂ）のパルスシーケンスを用いて取得される画像データである。ディフュージョンテンソル画像情報２００と同様に、キャビネット部２０の記憶部１６６から、画像処理装置４０の記憶部１７７に転送される。
【００５３】
テンソル行列算出手段２１０は、ディフュージョンテンソル画像情報２００およびＴ_２強調画像情報２０２から画素ごとの拡散テンソルを求める演算である。各座標軸ごとのＭＰＧ勾配の組み合わせ（ｂｘ、ｂｙ、ｂｚ）と、画素値Ｓと，拡散テンソルＤαβ（α、β∈ｘ、ｙ、ｚ）の間には、
【００５４】
【数１】

【００５５】
の関係がある。ここで、Ｓ_０は、Ｔ_２強調画像の画素値である。また、ｂは、ｂ値とよばれ、
ｂ＝γ^２δ^２（Δ―δ／３）
で定義される。ここで、γは、磁気回転比、δおよびΔは、図２（Ｂ）に示したＭＰＧ勾配の幅および間隔である。
【００５６】
拡散テンソルＤαβ（α、β∈ｘ、ｙ、ｚ）は、式（１）を用いて求められるが、テンソルの性質から独立要素は６つに限られるので、ＭＰＧ勾配を変化させた組み合わせ（ｂｘ、ｂｙ、ｂｚ）のディフュージョンテンソル画像情報２００は、最低６つ必要となる。また、それ以上のＭＰＧ勾配の組み合わせを用いた場合には、最小二乗法を用いて、拡散テンソルＤαβ（α、β∈ｘ、ｙ、ｚ）の値を高精度化する。
【００５７】
行列対角化手段２２０は、行列で表現される拡散テンソルＤαβを対角化し、固有ベクトル（ｖｅｃｔｏｒ）、固有値を求める演算である。図５にこの演算を例示した。ここで、拡散テンソルＤαβは、３行３列の行列であるので、３つの固有ベクトルと３つの固有値λ_ｉ（ｉ＝１，２，３）有する。この演算は、通常の行列の対角化により行われ，３つの固有ベクトルを要素として含む正則行列Ｑを求める（斉藤正彦著、線形代数入門、東京大学出版会、第５章参照）。
【００５８】
指標算出手段２３０は、拡散テンソルＤαβの固有値から、拡散テンソルの異方性の大きさを示す指標値を算出する。この指標としては、例えばフラクショナルアニソトロピー（ｆｒａｃｔｉｏｎａｌ　ａｎｉｓｏｔｒｏｐｙ：以下ＦＡ値と称する）が用いられる。ＦＡ値は、
【００５９】
【数２】

【００６０】
で定義される。ここで、λ_ｉは、拡散テンソルの３つの固有値の内の１つを現す。また、Ｉは、
Ｉ＝λ_１＋λ_２＋λ_３
で定義される。ＦＡ値によれば、３つの固有値λ_ｉが大きく異なれば、大きな値となり、３つの固有値λ_ｉが等しくなれば、零となる。すなわち、ＦＡ値が大きい場合は、異方性が大きく、ＦＡ値が小さい場合には、異方性が小さい。
【００６１】
閾値手段２４０は、操作部１９１から設定された閾値を用いて、画像データの中から、閾値を越える画素値の領域を抽出する手段である。この抽出領域は、記憶部１７７に保存され、画像処理等が施される。
【００６２】
領域拡張手段２５０は、画素ごとに有するパラメータ値、本実施の形態では最大固有ベクトルの向き、に基づいて、３×３×３画素の立方体状のマスク（ｍａｓｋ）領域を設定し、この中心画素の周辺に位置する画素領域の中から、中心画素の最大固有ベクトルの向きに所定範囲内で近似する画素を選択する。その後、選択された画素を新たな中心画素とするマスク領域を設定し、この周辺画素に対して同様の操作を繰り返し、中心画素の抽出を順次行う。ここで、新たな中心画素の抽出では、既に抽出された中心画素は排除されるので、順次抽出領域は拡張される。
【００６３】
図６に、平面上の画素に対して行われる領域拡張手段２５０の例を示した。図６（Ａ）には、３×３画素の正方形状のマスク領域と、この中心画素Ｐおよびその周辺画素の最大固有ベクトルが図示されている。なお、平面上の場合には、１つ中心画素に対して、８個の周辺画素が存在する。また、実際に行われる、立方体状のマスク領域を用いた領域拡張手段２５０の場合には、１つの中心画素に対して、２６個の周辺画素が存在する。
【００６４】
ここで、中止画素Ｐの最大固有ベクトルの向きと、所定範囲内で近似する新たな中心画素Ｐ′が選択される。そして、中心画素Ｐ′を中心にした新たなマスク領域が設定される（図６（Ｂ））。その後、図６（Ｃ）では、新たなマスク領域の周辺画素の中から、所定範囲内で近似する新たな中心画素Ｐ″が選択される。この様にして、中心画素Ｐ，Ｐ′、Ｐ″が順次抽出される。また、立方体状のマスク領域を用いる３次元的な領域拡張手段２５０でも全く同様である。
【００６５】
ボリュームレンダリング手段２６０は、抽出された立体的な繊維状の画像情報を、表示部１８１に平面的な画像として表示する際に、陰影処理、補間処理等の手法により、立体感を持たせ、オペレータの理解を助け、より見易くするものである。
【００６６】
つづいて、繊維状画像の抽出を行う画像処理装置４０の動作について図７を用いて説明する。まず、画像処理装置４０は、操作コンソール部３０から、Ｔ_２強調画像情報２０２およびこの画像のディフュージョンテンソル画像情報２００を取得する（ステップＳ６０１および６０２）。図９（Ａ）は、図２（Ａ）に示した頭部スライス面のＴ_２強調画像の表示画面の一例である。特に、神経繊維が、主として表示画面と平行して走査する脳梁部分および表示画面と垂直方向に走査する錐体路部分が描出されている。
【００６７】
そして、これらの画像情報から、テンソル行列算出手段２１０を用いて、画素ごとの拡散テンソル行列を算出する（ステップＳ６０３）。この拡散テンソル行列は、３行３列の行列からなり、各要素は、方向により異なる拡散係数の成分を現している。
【００６８】
その後、画素ごとに、拡散テンソル行列を、行列対角化手段２２０を用いて対角化し、その過程で、固有値、固有ベクトルを算出する（ステップＳ６０４）。この固有値、固有ベクトルの物理的意味を、図８に示す。すなわち、図８に示す神経繊維内の点においては、神経繊維の進行方向と、これと垂直をなす方向とで水分子の拡散係数が異なり、任意方向の拡散係数は、進行方向およびこれと垂直をなす拡散係数の線形結合により求めることができる。
【００６９】
この様な物理系においては、拡散テンソル行列は、対角化が可能であり、固有ベクトルの向きは、神経繊維の進行方向あるいはこれと垂直をなす方向を示している。また、神経繊維の進行方向の拡散係数Ｄ‖が、これと垂直をなす方向の拡散係数Ｄ⊥よりも大きくなることが知られている（前記阿部修等：ＭＲＩ：拡散画像；画像診断、Ｖｏｌ１９、Ｎｏ．１０、１９９９参照）。このことから、拡散係数である固有値が最大となる固有ベクトルの向きは、神経繊維の進行方向を指し示している。
【００７０】
なお、神経繊維が存在しない均一組織では、拡散係数は、すべての方向で等しいものとなり、固有値は、すべて等しくなる。
図７に戻り、ステップＳ６０４で求めた固有値から、指標算出手段２３０を用いて、ＦＡ値を算出する（ステップＳ６０５）。神経繊維が存在する領域では、拡散係数である固有値は大きく異なり、他方、神経繊維が存在しない均一な領域では、固有値は等しくなる。ここで、ＦＡ値は、前述した様に、固有値の異方性を示す指標である。この指標を用いることにより、神経繊維が存在する領域では、ＦＡ値は大きくなり、神経繊維が存在しない領域では、ＦＡ値は小さくなるので、神経繊維が存在する領域を識別することができる。
【００７１】
その後、ステップＳ６０５で算出された画素ごとのＦＡ値から、閾値手段２４０を用いて、神経繊維領域の抽出を行う（ステップＳ６０６）。ここでは、神経繊維領域において、ＦＡ値が高くなることを用いて、ＦＡ値に閾値を設け、この閾値を基準として、この基準を超える領域の抽出を行う。ここで、ＦＡ値の閾値としては、例えば０．６が設定され、ＦＡ値≧０．６の領域が抽出される。図９（Ｂ）に、図９（Ａ）のＴ_２強調画像から、神経繊維領域である脳梁および錐体路が、ＦＡ値をもとに抽出された様子を模式的に示した。なお、図では分かり易くするため平面的な抽出を示したが、本実施の形態では、ディフュージョンテンソル画像情報２００の撮像部位全体から、立体的な神経繊維領域が抽出される。
【００７２】
その後、ステップＳ６０６で抽出された神経繊維領域の中から、異方性の高い画素を始点として、領域拡張手段２５０を用いて、神経繊維の抽出を行う（ステップＳ６０７）。ここで、異方性の高い始点画素を、例えばＦＡ値≧０．８の領域からオペレータが選択する。そして、領域拡張手段２５０を用いる際に、中心画素を基点として周辺画素を選択するパラメータとしては、中心画素と周辺画素との結合指数（コネクティビティ：ｃｏｎｎｅｃｔｉｖｉｔｙとも称される）が用いられる。この結合指数には、中心画素と周辺画素の最大固有ベクトルのなす方向余弦の角度が用いられる。
【００７３】
この中心画素と周辺画素の最大固有ベクトルのなす角度は、中心画素と周辺画素の最大固有ベクトルの内積から求められる。中心画素の最大固有ベクトルをν、周辺画素の最大固有ベクトルをν′とすると、中心画素と周辺画素の最大固有ベクトルのなす角度θは、内積と
ν・ν′＝｜ν｜｜ν′｜ｃｏｓθ
の関係式から求められる。なお、ステップＳ６０４で算出される固有ベクトルが規格化されている場合には、｜ν｜＝｜ν′｜＝１となるので、ベクトルνおよびν′の内積は、ベクトルνおよびν′の方向余弦ｃｏｓθとなる。そして、この角度θに閾値を設け、閾値以下の周辺画素、すなわち最大固有ベクトルが概ね同じ方向を有する周辺画素を新たな中心画素とする。これにより、神経繊維の進行方向の画素の抽出を、順次行う。ここで、この角度θの閾値としては、例えばν・ν′≧０．９９６、すなわちθ≦５．１３度が経験上設定される。
【００７４】
図９（Ｃ）に、図９（Ｂ）の神経繊維領域から抽出された２次元的な神経繊維画像の例を模式的に示した。図９（Ｂ）に示された神経繊維領域である錐体路および脳梁は、錐体路では概ね表示画面に垂直方向に神経繊維が走行し、脳梁では概ね表示画面に沿って神経繊維が走行している。このため、主として、脳梁部分の神経繊維が、表示画面上の領域拡張手段２５０により抽出される。
【００７５】
その後、ステップＳ６０７で抽出した３次元的な神経繊維画像情報を、ボリュームレンダリング手段２６０を用いて、平面上に投影するボリュームレンダリングを行う（ステップＳ６０８）。これにより、表示部１８１のＣＲＴ画面に表示する際に、オペレータにとって見易く、しかも神経繊維の構造が理解しやすいものとなる。そして、この結果を、表示部１８１に表示して（ステップＳ６０９）、本処理を終了する。
【００７６】
上述してきたように、本実施の形態では、ディフュージョンテンソル画像情報２００およびＴ_２強調画像情報２０２から、拡散テンソル行列を算出し、この拡散テンソル行列の固有値および固有ベクトルを行列の対角化により求め、この固有値から拡散係数の異方性の指標であるＦＡ値を算出し、このＦＡ値から神経繊維領域を抽出し、さらに神経繊維の進行方向である最大固有ベクトルの方向から、領域拡張手段２５０を用いて、神経繊維領域から順次神経繊維画素を抽出することとしているので、短時間で、しかも個人差もなく神経繊維画像情報の抽出を行うことができる。
【００７７】
また、本実施の形態では、操作コンソール部３０から取得したＴ_２強調画像情報２０２およびこの画像のディフュージョンテンソル画像情報２００に対して、拡散テンソル行列の算出を行ったが、フィルター（ｆｉｌｔｅｒ）処理によるスムージング（ｓｍｏｏｔｈｉｎｇ）をこれら画像情報に施した後に拡散テンソル行列の算出を行うこともできる。
【００７８】
【発明の効果】
以上説明したように、本発明によれば、イメージングパルスシーケンスを実行して被検体の同一撮像部位のディフュージョンテンソル画像情報およびＴ_２強調画像情報を取得し、このディフュージョンテンソル画像情報およびＴ_２強調画像情報から画素ごとの拡散テンソル行列を算出し、この拡散テンソル行列の固有値および固有ベクトルを算出し、この固有値から画素位置における拡散係数の異方性の指標値を算出し、この指標値に基づいて、撮像部位から異方性領域を抽出し、この異方性領域内の指標値が閾値を越える領域に始点画素を選定し、画素ごとの固有ベクトルの中から、最大の固有値を有する最大固有ベクトルを選択し、始点画素から、最大固有ベクトルのベクトル情報を用い、繊維状画像を順次抽出することとしているので、神経繊維の進行方向を示す最大固有ベクトルを算出し、このベクトルの方向に沿った神経繊維画素を順次，自動で抽出することができ、ひいては、神経繊維画素の抽出に係る時間を大幅に削減し、さらに、抽出における個人差をなくすことができる。
【図面の簡単な説明】
【図１】磁気共鳴撮像システムの全体構成を示すブロック図である。
【図２】実施の形態の撮像部位およびパルスシーケンスを示す図である。
【図３】実施の形態のディフュージョンテンソル画像情報を示す図である。
【図４】実施の形態の画像処理装置の機能ブロック図である。
【図５】拡散テンソル行列の対角化を示す図である。
【図６】領域拡張手段の動作を示す図である。
【図７】実施の形態の画像処理装置の動作を示すフローチャートである。
【図８】神経繊維内の拡散係数の大きさおよび向きを示す図である。
【図９】神経繊維画像の抽出を模式的に示した図である。
【符号の説明】
１　被検体
２０　キャビネット部
３０　操作コンソール部
４０　画像処理装置
１００　マグネット部
１０２　主磁場コイル部
１０６　勾配コイル部
１０８　ＲＦコイル部
１２０　クレードル
１３０　勾配駆動部
１４０　ＲＦ駆動部
１５０　データ収集部
１５５　データ処理部
１６０　スキャンコントローラ部
１７０　制御部
１７１　画像処理部
１６６，１７６、１７７　記憶部
１８０、１８１　表示部
１９０、１９１　操作部
２００　ディフュージョンテンソル画像情報
２０２　Ｔ_２強調画像情報
２１０　テンソル行列算出手段
２２０　行列対角化手段
２３０　指標算出手段
２４０　閾値手段
２５０　領域拡張手段
２６０　ボリュームレンダリング手段[0001]
TECHNICAL FIELD OF THE INVENTION
The present invention relates to a fibrous image extraction method, an image processing apparatus, and a magnetic resonance imaging system, which are extracted from an imaging site using diffusion tensor image information.
[0002]
[Prior art]
In recent years, in the field of magnetic resonance imaging apparatuses, in addition to conventional morphological image information of a subject, a functional image obtained by imaging the function of the subject has been clinically used. Among the functional images, the diffusion-weighted image is a method of calculating a diffusion coefficient in a tissue from a diffusion tensor image and forming an image. It is known that the anisotropy of the diffusion coefficient indicates the running direction of the white matter nerve fibers in the brain or spinal cord (Osamu Abe et al .: MRI: diffusion image; diagnostic imaging, Vol19, No. 10, 1999).
[0003]
Therefore, a region where a nerve fiber exists is selected based on the magnitude of the anisotropy of the diffusion coefficient, and a nerve fiber image indicating the running direction of the nerve fiber is extracted from the selected region.
[0004]
[Problems to be solved by the invention]
However, according to the above-described conventional technology, a great deal of time and specialized knowledge are required to extract only the nerve fiber image. In other words, the operation of extracting the nerve fiber image from the selected region requires the doctor to manually select and select each pixel based on anatomical knowledge.
[0005]
In particular, when a doctor manually extracts only nerve fiber images, individual differences depending on anatomical knowledge occur, and it takes a lot of time to select and remove a large amount of image information, This is a cause of human error caused by this long operation.
[0006]
From these facts, it is important how to realize a fibrous image extraction method, an image processing apparatus, and a magnetic resonance imaging system that can automatically and quickly extract only a nerve fiber image.
[0007]
SUMMARY OF THE INVENTION The present invention has been made to solve the above-described problem of the related art, and provides a fibrous image extraction method, an image processing apparatus, and a magnetic resonance imaging system capable of automatically and quickly extracting only a nerve fiber image. The purpose is to do.
[0008]
[Means for Solving the Problems]
In order to solve the above-mentioned problems and achieve the object, a fibrous image extraction method according to a first aspect of the present invention executes an imaging pulse sequence to execute diffusion tensor image information and T₂Acquiring the enhanced image information, the diffusion tensor image information and the T₂A diffusion tensor matrix for each pixel is calculated from the emphasized image information, an eigenvalue and an eigenvector of the diffusion tensor matrix are calculated, and an index value of anisotropy of a diffusion coefficient at the pixel position is calculated from the eigenvalue. Based on the extracted anisotropic region from the imaging region, the index value in the anisotropic region is selected as a starting pixel in a region exceeding a threshold, from among the eigenvectors for each pixel, the maximum of the A maximum eigenvector having an eigenvalue is selected, and a fibrous image is sequentially extracted from the starting pixel using vector information of the maximum eigenvector.
[0009]
According to the first aspect of the present invention, the diffusion tensor image information and the T₂The enhanced tensor image information is acquired, and the diffusion tensor image information and T₂A diffusion tensor matrix for each pixel is calculated from the emphasized image information, an eigenvalue and an eigenvector of the diffusion tensor matrix are calculated, an index value of anisotropy of a diffusion coefficient at a pixel position is calculated from the eigenvalue, and based on the index value, Then, an anisotropic region is extracted from the imaging region, a starting pixel is selected in a region where the index value in the anisotropic region exceeds a threshold, and a maximum eigenvector having a maximum eigenvalue is determined from eigenvectors for each pixel. From the start pixel, the vector image of the maximum eigenvector is used to sequentially extract fibrous images from the start pixel.Therefore, the maximum eigenvector indicating the traveling direction of the nerve fiber is calculated, and the nerve fiber pixels along the direction of this vector are calculated. Can be sequentially and automatically extracted.
[0010]
Further, the image processing apparatus according to the invention of the second aspect provides the diffusion tensor image information and the T₂Storage means for storing enhanced image information, the diffusion tensor image information and the T₂Tensor matrix calculating means for calculating a diffusion tensor matrix for each pixel from the emphasized image information; matrix diagonalization means for calculating eigenvalues and eigenvectors of the diffusion tensor matrix; and anisotropy of a diffusion coefficient at the pixel position from the eigenvalues An index processing unit that calculates an index value of, and a threshold unit that extracts a high anisotropic region from the imaging region based on the index value, an image processing apparatus including: A starting pixel is selected in an area where the index value exceeds a threshold, and among the eigenvectors for each pixel, a maximum eigenvector having the largest eigenvalue is selected.From the starting pixel, vector information of the maximum eigenvector is used. And a region expanding means for extracting a fibrous image.
[0011]
According to the second aspect of the present invention, the starting pixel is selected in a region where the index value in the anisotropic region exceeds the threshold value by using the region expanding means, and the maximum eigenvalue is determined from the eigenvectors for each pixel. Since the maximum eigenvector is selected and the fibrous image is extracted from the starting pixel using the vector information of the maximum eigenvector, the maximum eigenvector indicating the traveling direction of the nerve fiber is calculated, and the pixel in that direction is calculated as the area. The nerve fiber image can be automatically extracted by performing the extraction by the extension means.
[0012]
Further, an image processing apparatus according to a third aspect of the invention is characterized in that a fractional anisotropy value calculated from the eigenvalue is used as the index value.
[0013]
According to the third aspect of the invention, since the fractional anisotropy value calculated from the eigenvalue is used as the index value, the anisotropy of the diffusion coefficient can be quantitatively evaluated.
[0014]
An image processing apparatus according to a fourth aspect of the present invention is characterized in that the threshold value means extracts a pixel region in which the fractional anisotropy value exceeds a threshold value.
[0015]
According to the fourth aspect of the invention, since the threshold value means extracts a pixel region having a fractional anisotropy value exceeding the threshold value, a nerve fiber region having high anisotropy is extracted by the fractional anisotropy value. be able to.
[0016]
An image processing apparatus according to a fifth aspect of the invention is characterized in that the tensor matrix calculating means is optimized using a least squares method.
According to the fifth aspect of the invention, the tensor matrix calculation means is optimized using the least squares method, so that a highly accurate diffusion tensor matrix can be obtained.
[0017]
The image processing apparatus according to a sixth aspect of the present invention is the image processing apparatus, wherein the area expansion unit calculates a connection index indicating a connection between a central pixel and a peripheral pixel of a mask area for each of the peripheral pixels from the vector information. A new mask area is set with the peripheral pixel whose coupling index does not exceed a threshold value as a new center pixel, and area expansion is performed sequentially.
[0018]
According to the invention of the sixth aspect, the region expanding means calculates a connection index indicating a connection between the central pixel and the peripheral pixels of the mask region for each peripheral pixel from the vector information, and the connection index exceeds the threshold value. Since a new mask area is set with no peripheral pixel as a new central pixel, and the area is expanded sequentially, a peripheral pixel made of the same nerve fiber as that of the central pixel is selected, and area expansion is performed. Thus, the whole nerve fiber can be extracted.
[0019]
An image processing apparatus according to a seventh aspect of the invention is characterized in that an angle of a direction cosine calculated from an inner product of the vector information of the central pixel and the peripheral pixels is used as the connection index.
[0020]
According to the seventh aspect of the invention, since the angle of the direction cosine calculated from the inner product of the vector information of the center pixel and the peripheral pixels is used as the joint index, pixels having the same direction of the maximum eigenvector are determined. , Binding index.
[0021]
A magnetic resonance imaging system according to an eighth aspect of the present invention forms a static magnetic field, forms a gradient magnetic field, transmits or receives an RF signal, controls the gradient magnetic field and the RF signal, Tensor image information and T₂Magnetic resonance imaging apparatus for acquiring enhanced image information, and the diffusion tensor image information and T of the magnetic resonance imaging apparatus₂Transfer means for transferring the emphasized image information; diffusion tensor image information transferred from the transfer means;₂A magnetic resonance imaging system comprising: an image processing device that performs image processing on the emphasized image information, wherein the image processing device includes the diffusion tensor image information and the T₂Tensor matrix calculating means for calculating a diffusion tensor matrix for each pixel from the emphasized image information; matrix diagonalization means for calculating eigenvalues and eigenvectors of the diffusion tensor matrix; and anisotropy of a diffusion coefficient at the pixel position from the eigenvalues Index calculating means for calculating an index value of, a threshold means for extracting an anisotropic region from the imaging region based on the index value, and a starting point in a region in the anisotropic region where the index value exceeds a threshold value Area expansion means for selecting a pixel, selecting the largest eigenvector having the largest eigenvalue from among the eigenvectors for each pixel, and extracting a fibrous image from the starting pixel using the vector information of the largest eigenvector. And the following.
[0022]
According to the eighth aspect of the present invention, in the image processing apparatus, the tensor matrix calculating means uses the diffusion tensor image information and T₂A diffusion tensor matrix for each pixel is calculated from the emphasized image information, an eigenvalue and an eigenvector of the diffusion tensor matrix are calculated by matrix diagonalization means, and an anisotropy of a diffusion coefficient at a pixel position is calculated from the eigenvalue by an index calculation means. The anisotropic region is extracted from the imaging region by the threshold means based on the index value, and the starting point is set to the region where the index value in the anisotropic region exceeds the threshold by the region expanding means. A pixel is selected, a maximum eigenvector having a maximum eigenvalue is selected from eigenvectors for each pixel, and a fibrous image is extracted from a starting pixel using vector information of the maximum eigenvector. Is calculated by calculating the maximum eigenvector indicating the traveling direction of the nerve fiber, and extracting the pixels in that direction by the region expanding means, thereby automatically extracting the nerve fiber image. It can be.
[0023]
BEST MODE FOR CARRYING OUT THE INVENTION
Exemplary embodiments of a fibrous image extraction method, an image processing apparatus, and a magnetic resonance imaging system according to the present invention will be described below with reference to the accompanying drawings. Note that the present invention is not limited to this.
[0024]
First, the overall configuration of the magnetic resonance imaging system according to the present embodiment will be described. FIG. 1 shows a block diagram of the magnetic resonance imaging system. As shown in FIG. 1, the present system includes a magnet unit 100, a cabinet unit 20, an operation console unit 30, and an image processing device 40.
[0025]
The magnet unit 100 has a main magnetic field coil unit 102, a gradient coil unit 106, and an RF (radio frequency) coil unit 108. Each of these coil portions has a substantially cylindrical shape and is arranged coaxially with each other. The subject 1 to be imaged in a recumbent state on a cradle 120 is arranged at the center of a substantially cylindrical imaging space (bore) of the magnet unit 100.
[0026]
The main magnetic field coil unit 102 forms a static magnetic field in the internal space of the magnet unit 100. The direction of the static magnetic field is substantially parallel to the axial direction of the cylindrical coil arranged coaxially. That is, a so-called horizontal magnetic field is formed. The main magnetic field coil unit 102 is configured using, for example, a superconducting coil. In addition, you may comprise using not only a superconducting coil but a normal conduction coil.
[0027]
The gradient coil unit 106 includes three independent gradient coils (not shown), that is, three gradient coils for giving a linear gradient to the static magnetic field strength in the x, y, and z coordinate axis directions. Here, a slice axis, a phase encode axis, and a frequency encode axis are arbitrary coordinate axes that determine an imaging section and that are orthogonal to each other with respect to the x, y, and z coordinate axes while maintaining perpendicularity between them. It is also possible to have a slope.
[0028]
The gradient magnetic field in the slice axis direction is also called a slice gradient magnetic field, the gradient magnetic field in the phase encode axis direction is also called a phase encode gradient magnetic field, and the gradient magnetic field in the frequency encode axis direction is also called a frequency encode gradient magnetic field. Hereinafter, the gradient magnetic field is also simply referred to as a gradient.
[0029]
The RF coil unit 108 forms a high-frequency magnetic field for exciting a spin in the body of the subject 1 in the static magnetic field space. Hereinafter, forming a high-frequency magnetic field is also referred to as transmitting an RF excitation signal. The RF excitation signal is also called an RF pulse. The RF coil unit 108 receives the electromagnetic wave generated by the excited spin, that is, the magnetic resonance signal, together with the transmission of the RF excitation signal.
[0030]
The RF coil unit 108 has a transmitting coil and a receiving coil. The same coil may be used for the transmitting coil and the receiving coil, or dedicated coils may be used.
[0031]
The cabinet 20 includes a scan controller 160, a gradient driver 130, an RF driver 140, a data collector 150, a data processor 155, and a storage 166.
[0032]
The scan controller 160 controls the gradient driving unit 130, the RF driving unit 140, the data collection unit 150, and the data processing unit 155.
The gradient driving unit 130 is connected to the gradient coil unit 106. The gradient driving section 130 supplies a driving signal to the gradient coil section 106 to generate a gradient magnetic field. The gradient drive unit 130 has three drive circuits (not shown) corresponding to the three gradient coils in the gradient coil unit 106.
[0033]
The RF coil unit 108 is connected to the RF driving unit 140. The RF driving unit 140 supplies a driving signal to the RF coil unit 108 to transmit an RF pulse to excite spins in the subject 1.
[0034]
The data collection unit 150 is also connected to the RF coil unit 108. The data collection unit 150 captures the received signal received by the RF coil unit 108 by sampling, and collects it as digital data (digital @ data). The output side of the data collection unit 150 is connected to a storage unit 166, and these digital data are stored and stored in the storage unit 166 as raw data.
[0035]
The scan controller 160 is configured using, for example, a computer. The scan controller 160 has a memory (not shown). The memory stores a pulse sequence (pulse @ sequence), which is a control program (program) for the scan controller unit 160, and various data. The function of the scan controller 160 is realized by the computer executing a pulse sequence stored in the memory.
[0036]
A digitized received signal is stored as raw data in the storage unit 166, and the received signal is transferred to the data processing unit 155 and image reconstruction is performed upon completion of imaging. The data processing unit 155 is configured using a computer or the like, and performs image reconstruction by two-dimensional Fourier transform on raw data in K space, which is a magnetic resonance signal collected by the data collection unit 150, to generate image data. . Then, the image data is transferred to the storage unit 166 and stored.
[0037]
The control section 170 of the operation console section 30 is above the scan controller section 160 and the data processing section 155 of the cabinet section 20 and controls them. In addition, these operations are realized by the control unit 170 executing a program stored in the storage unit 176.
[0038]
The operation console unit 30 includes a control unit 170, a storage unit 176, a display unit 180, and an operation unit 190. The storage unit 176 is connected to the storage unit 166 of the cabinet unit 20 by a communication unit (not shown), and acquires the image data of the storage unit 166. The display unit 180 is configured with a graphic display (graphic display) or the like. The operation unit 190 includes a keyboard having a pointing device.
[0039]
The display unit 180 displays the image data of the storage unit 176 and various information. The operation unit 190 is operated by an operator and inputs various commands and information to the control unit 170. The operator operates the present apparatus interactively through the display unit 180 and the operation unit 190.
[0040]
The image processing device 40 includes a storage unit 177, an image processing unit 171, a display unit 181, and an operation unit 191. The storage unit 177 is connected to the storage unit 176 of the operation console unit 30 by communication means (not shown), and acquires image data of the storage unit 176. The image processing unit 171 performs image processing on the image data in the storage unit 177, and displays the image data on the display unit 181. The display unit 181 is configured by a graphic display or the like. The operation unit 191 is configured by a keyboard or the like having a pointing device.
[0041]
The display unit 181 displays the image data, the image processing data, and various information of the storage unit 177. The operation unit 191 is operated by an operator, and inputs various commands, information, and the like to the image processing unit 171. The operator interactively operates the apparatus through the display unit 181 and the operation unit 191.
[0042]
Here, before describing the functional block diagram and operation of the image processing apparatus 40, acquisition of a diffusion tensor image using an echo planar method will be described with reference to FIGS. First, the diffusion tensor image is image information of data serving as a basis for calculating a distribution of a diffusion coefficient of water protons (proton) existing in the subject 1. Here, the diffusion coefficient in the subject 1 has a different value depending on the direction even at the same pixel position, and is described by a tensor amount using a matrix.
[0043]
The operator sets the imaging surface of the subject 1 and the number of slices for imaging from the operation unit 190. FIG. 2A is a diagram showing an example in which the head of the subject 1 is sliced on N tilted axial imaging surfaces from a sagittal cross section. The slice plane and the number of slices are determined in consideration of the range of the nerve fiber region to be imaged and the direction of the nerve fiber.
[0044]
Thereafter, image information for each slice is acquired by using an echo planar method using a spin echo to which a motion detection gradient magnetic field (hereinafter, referred to as MPG gradient: motion {probing} gradient) is added. FIG. 2B shows an example of the pulse sequence of the echo planar method. FIG. 2B is a time chart (time @ chart) showing operations in three axis directions of a slice axis, a phase encode axis, and a frequency encode axis, which have a common time axis on the horizontal axis.
[0045]
The RF pulse in FIG. 2B excites a target slice plane of the subject 1 specified by the slice gradient by a 90-degree pulse, and further inverts the phase of the excited proton by a 180-degree pulse to produce a spin echo. Take out.
[0046]
The phase encoding gradient applies a phase change to the excited protons in the phase encoding axis direction, and extracts frequency component information corresponding to the phase change. Further, the frequency encoding gradient is also called a readout gradient, and a gradient is added in the frequency encoding axis direction, and the position inside the subject 1 is specified from the frequency component of the spin echo observed during the application of the gradient. When a slice plane that forms an angle with the body axis direction of the subject 1 as shown in FIG. 2A is selected, the phase encode gradient, the frequency encode gradient, and the slice gradient are x, y, It does not match the direction of the z-coordinate axis, and is distributed to the x, y, and z-coordinate axes by the same transformation as the vector coordinate transformation while maintaining the perpendicularity between them.
[0047]
The MPG gradient shown in FIG. 2B is a gradient having the same magnitude and the same time interval, which is disposed at a symmetric position with a 180-degree pulse interposed therebetween. Although the MPG gradient is shown in FIG. 2B on a different axis from the slice axis, the phase encode axis, and the frequency encode axis, in an actual scan, at least the slice axis, the phase encode axis, and the frequency encode axis are shown. The size is changed to one and added. As a result, the echo signal of water protons that have a random motion and do not have regularity and that contribute to the diffusion coefficient is reduced due to phase diffusion. Further, by changing the magnitude of the MPG gradient and the added axial direction, the anisotropy information of the diffusion coefficient is also obtained.
[0048]
Here, the magnitude of the MPG gradient applied in the x coordinate axis direction is bx, the magnitude of the MPG gradient applied in the y coordinate axis direction is by, and the magnitude of the MPG gradient applied in the z coordinate axis direction is bz. This gradient is arranged on the same time axis as the slice axis, the phase encode axis, the frequency encode axis, or the x, y, and z coordinate axes. Then, for each combination (bx, by, bz) of the magnitude of the MPG gradient of each of the x, y, and z coordinate axes, raw data of N slice planes is acquired.
[0049]
FIG. 3 schematically shows an example of diffusion tensor image information 200 obtained from imaging of m combinations of different MPG gradients. 2B is image data obtained by reconstructing an image of the N pieces of raw data, which is the number of slice planes, acquired by the pulse sequence of FIG. 2B for each combination of MPG gradients. In order to calculate a diffusion tensor matrix to be described later, it is necessary that m> 6. Also, combinations of different MPG gradients are indicated by a subscript m.
[0050]
Subsequently, a functional configuration of the image processing apparatus 40 will be described. FIG. 4 is a functional block diagram of the image processing apparatus 40. The image processing device 40 includes diffusion tensor image information 200,₂Enhanced image information 202, a tensor matrix calculation unit 210 as a calculation unit, a matrix diagonalization unit 220, an index calculation unit 230, a threshold unit 240 as a region setting unit, a region expansion unit 250, and a volume rendering (volume) as a display unit rendering means 260 and the like.
[0051]
The diffusion tensor image information 200 is transferred from the storage unit 166 of the cabinet unit 20 to the storage unit 176 of the operation console unit 30, and further transferred to the storage unit 177 of the image processing device 40. As shown in FIG. 3, the diffusion tensor image information 200 is composed of the number of slices, for example, N pieces of image data for each number of combinations of the magnitude of the MPG gradient, for example, m.
[0052]
T₂The emphasized image information 202 is image data obtained by using the pulse sequence of FIG. 2B when the magnitude of the MPG gradient of the diffusion tensor image information 200 becomes zero on all coordinate axes. Similarly to the diffusion tensor image information 200, the information is transferred from the storage unit 166 of the cabinet unit 20 to the storage unit 177 of the image processing device 40.
[0053]
The tensor matrix calculating means 210 calculates the diffusion tensor image information 200 and T₂This is an operation for obtaining a diffusion tensor for each pixel from the emphasized image information 202. Between the combination (bx, by, bz) of the MPG gradient for each coordinate axis, the pixel value S, and the diffusion tensor Dαβ (α, β∈x, y, z),
[0054]
(Equation 1)

[0055]
There is a relationship. Where S₀Is T₂This is the pixel value of the emphasized image. B is called a b value,
b = γ²δ²(Δ-δ / 3)
Is defined by Here, γ is the gyromagnetic ratio, and δ and Δ are the width and interval of the MPG gradient shown in FIG.
[0056]
The diffusion tensor Dαβ (α, β∈x, y, z) is obtained by using the equation (1). However, since the number of independent elements is limited to six due to the nature of the tensor, the combination (bx , By, bz) of the diffusion tensor image information 200 is required at least six. When a combination of more MPG gradients is used, the value of the diffusion tensor Dαβ (α, β∈x, y, z) is made highly accurate using the least squares method.
[0057]
The matrix diagonalizing means 220 is an operation for diagonalizing the diffusion tensor Dαβ expressed by the matrix to obtain an eigenvector (vector) and an eigenvalue. FIG. 5 illustrates this calculation. Here, since the diffusion tensor Dαβ is a matrix of three rows and three columns, three eigenvectors and three eigenvalues λ_i(I = 1, 2, 3). This operation is performed by ordinary diagonalization of the matrix to obtain a regular matrix Q including three eigenvectors as elements (see Masahiko Saito, Introduction to Linear Algebra, University of Tokyo Press, Chapter 5).
[0058]
The index calculator 230 calculates an index value indicating the magnitude of the anisotropy of the diffusion tensor from the eigenvalue of the diffusion tensor Dαβ. As this index, for example, fractional anisotropy (hereinafter, referred to as FA value) is used. FA value is
[0059]
(Equation 2)

[0060]
Is defined by Where λ_iRepresents one of the three eigenvalues of the diffusion tensor. Also, I
I = λ₁+ Λ₂+ Λ₃
Is defined by According to the FA value, three eigenvalues λ_iAre greatly different, three eigenvalues λ are obtained._iIf are equal, it is zero. That is, when the FA value is large, the anisotropy is large, and when the FA value is small, the anisotropy is small.
[0061]
The threshold unit 240 is a unit that extracts a region having a pixel value exceeding the threshold from the image data using the threshold set from the operation unit 191. This extracted area is stored in the storage unit 177 and subjected to image processing and the like.
[0062]
The region expanding means 250 sets a 3 × 3 × 3 pixel cubic mask region based on the parameter value of each pixel, in this embodiment, the direction of the largest eigenvector, and From the pixel regions located in the periphery, a pixel that approximates the direction of the maximum eigenvector of the central pixel within a predetermined range is selected. After that, a mask area is set with the selected pixel as a new central pixel, and the same operation is repeated for the peripheral pixels to sequentially extract the central pixel. Here, in the extraction of a new center pixel, the already extracted center pixel is excluded, so that the extraction region is sequentially expanded.
[0063]
FIG. 6 shows an example of the area expanding means 250 performed on pixels on a plane. FIG. 6A shows a square mask area of 3 × 3 pixels, and the maximum eigenvector of the central pixel P and its surrounding pixels. In the case of a plane, there are eight peripheral pixels for one central pixel. Also, in the case of the area expanding means 250 using a cubic mask area, which is actually performed, there are 26 peripheral pixels for one central pixel.
[0064]
Here, a new center pixel P ′ that is similar to the direction of the largest eigenvector of the stop pixel P within a predetermined range is selected. Then, a new mask area centering on the center pixel P 'is set (FIG. 6B). Thereafter, in FIG. 6C, a new central pixel P ″ that is approximated within a predetermined range is selected from the peripheral pixels of the new mask area. In this manner, the central pixels P, P ′, P ″ Are sequentially extracted. The same applies to the three-dimensional region expanding means 250 using a cubic mask region.
[0065]
When displaying the extracted three-dimensional fibrous image information as a two-dimensional image on the display unit 181, the volume rendering means 260 gives a three-dimensional effect by a method such as shading processing and interpolation processing. It will help you understand and make it easier to see.
[0066]
Next, the operation of the image processing apparatus 40 for extracting a fibrous image will be described with reference to FIG. First, the image processing device 40 sends the T₂The emphasis image information 202 and the diffusion tensor image information 200 of this image are acquired (steps S601 and 602). FIG. 9 (A) shows the T of the head slice plane shown in FIG. 2 (A).₂It is an example of a display screen of an emphasized image. In particular, a corpus callosum portion in which nerve fibers scan mainly in parallel with the display screen and a pyramidal tract portion in which the nerve fiber scans in a direction perpendicular to the display screen are depicted.
[0067]
Then, a diffusion tensor matrix for each pixel is calculated from the image information by using the tensor matrix calculation unit 210 (step S603). This diffusion tensor matrix is composed of a matrix of 3 rows and 3 columns, and each element represents a component of a diffusion coefficient that differs depending on the direction.
[0068]
Then, for each pixel, the diffusion tensor matrix is diagonalized using the matrix diagonalization means 220, and in the process, eigenvalues and eigenvectors are calculated (step S604). FIG. 8 shows the physical meaning of these eigenvalues and eigenvectors. That is, at the points in the nerve fiber shown in FIG. 8, the diffusion coefficient of the water molecule is different between the traveling direction of the nerve fiber and the direction perpendicular thereto, and the diffusion coefficient in an arbitrary direction is different from the traveling direction and the vertical direction. Can be obtained by a linear combination of diffusion coefficients.
[0069]
In such a physical system, the diffusion tensor matrix can be diagonalized, and the direction of the eigenvector indicates the traveling direction of the nerve fiber or the direction perpendicular thereto. It is also known that the diffusion coefficient D‖ in the traveling direction of the nerve fiber is larger than the diffusion coefficient D⊥ in the direction perpendicular thereto (Osamu Abe et al .: MRI: diffusion image; image diagnosis, Vol19) , No. 10, 1999). From this, the direction of the eigenvector at which the eigenvalue that is the diffusion coefficient is the maximum indicates the traveling direction of the nerve fiber.
[0070]
In a uniform tissue where no nerve fiber exists, the diffusion coefficient is equal in all directions, and the eigenvalues are all equal.
Returning to FIG. 7, an FA value is calculated from the eigenvalue obtained in step S604 by using the index calculating means 230 (step S605). In the region where the nerve fiber exists, the eigenvalue which is the diffusion coefficient is greatly different, while in the uniform region where the nerve fiber does not exist, the eigenvalue becomes equal. Here, the FA value is an index indicating the anisotropy of the eigenvalue, as described above. By using this index, the FA value increases in an area where nerve fibers exist, and the FA value decreases in an area where nerve fibers do not exist, so that an area where nerve fibers exist can be identified.
[0071]
Thereafter, a nerve fiber region is extracted from the FA value for each pixel calculated in step S605 by using the threshold means 240 (step S606). Here, in the nerve fiber region, a threshold value is provided for the FA value by using the fact that the FA value becomes high, and an area exceeding this reference is extracted based on the threshold value. Here, for example, 0.6 is set as the threshold value of the FA value, and an area where FA value ≧ 0.6 is extracted. FIG. 9B shows a graph of T in FIG.₂From the enhanced image, a manner in which the corpus callosum and the pyramidal tract, which are nerve fiber regions, were extracted based on the FA value is schematically shown. Although a planar extraction is shown in the figure for simplicity, in the present embodiment, a three-dimensional nerve fiber region is extracted from the entire imaging region of the diffusion tensor image information 200.
[0072]
After that, nerve fibers are extracted from the nerve fiber region extracted in step S606, using the region having high anisotropy as a starting point, using the region expanding means 250 (step S607). Here, the operator selects a highly anisotropic starting pixel from, for example, an area where the FA value is ≥0.8. When using the area expansion unit 250, a coupling index (also referred to as connectivity) between the central pixel and the peripheral pixels is used as a parameter for selecting peripheral pixels based on the central pixel. The angle of the cosine of the direction formed by the maximum eigenvector of the central pixel and the peripheral pixels is used as the coupling index.
[0073]
The angle formed by the maximum eigenvector of the central pixel and the peripheral pixel is obtained from the inner product of the maximum eigenvector of the central pixel and the peripheral pixel. Assuming that the maximum eigenvector of the central pixel is ν and the maximum eigenvector of the peripheral pixels is ν ′, the angle θ between the central pixel and the maximum eigenvector of the peripheral pixels is
ν · ν '= | ν || ν' | cos θ
Is obtained from the relational expression. If the eigenvectors calculated in step S604 are normalized, | ν | = | ν ′ | = 1, so the inner product of the vectors ν and ν ′ is the direction cosine of the vectors ν and ν ′. cos θ. Then, a threshold value is provided for this angle θ, and peripheral pixels that are equal to or smaller than the threshold value, that is, peripheral pixels whose maximum eigenvectors have substantially the same direction, are set as new central pixels. Thus, the extraction of the pixels in the traveling direction of the nerve fiber is sequentially performed. Here, as a threshold value of the angle θ, for example, ν · ν ′ ≧ 0.996, that is, θ ≦ 5.13 degrees is empirically set.
[0074]
FIG. 9C schematically shows an example of a two-dimensional nerve fiber image extracted from the nerve fiber region in FIG. 9B. In the pyramidal tract and corpus callosum which are the nerve fiber regions shown in FIG. 9 (B), the nerve fiber runs in a direction substantially perpendicular to the display screen in the pyramidal tract, and the nerve fibers generally follow the display screen in the corpus callosum. Is running. For this reason, the nerve fibers in the corpus callosum are mainly extracted by the region expanding means 250 on the display screen.
[0075]
Thereafter, volume rendering is performed to project the three-dimensional nerve fiber image information extracted in step S607 on a plane using the volume rendering means 260 (step S608). Thus, when displayed on the CRT screen of the display unit 181, the operator can easily see and easily understand the structure of nerve fibers. Then, the result is displayed on the display unit 181 (step S609), and the process ends.
[0076]
As described above, in the present embodiment, diffusion tensor image information 200 and T₂A diffusion tensor matrix is calculated from the emphasized image information 202, an eigenvalue and an eigenvector of the diffusion tensor matrix are obtained by diagonalization of the matrix, and an FA value which is an index of anisotropy of a diffusion coefficient is calculated from the eigenvalue. Since the nerve fiber region is extracted from the FA value, and the nerve fiber pixel is sequentially extracted from the nerve fiber region using the region expanding means 250 from the direction of the maximum eigenvector which is the direction of the nerve fiber, In addition, the nerve fiber image information can be extracted without individual differences.
[0077]
Further, in the present embodiment, the T₂The diffusion tensor matrix was calculated for the emphasized image information 202 and the diffusion tensor image information 200 of this image. The diffusion tensor matrix was calculated after applying smoothing (smoothing) to these image information by filter processing. Can also be performed.
[0078]
【The invention's effect】
As described above, according to the present invention, the imaging pulse sequence is executed to execute diffusion tensor image information and T₂The enhanced tensor image information is acquired, and the diffusion tensor image information and T₂A diffusion tensor matrix for each pixel is calculated from the emphasized image information, an eigenvalue and an eigenvector of the diffusion tensor matrix are calculated, an index value of anisotropy of a diffusion coefficient at a pixel position is calculated from the eigenvalue, and based on the index value, Then, an anisotropic region is extracted from the imaging region, a starting pixel is selected in a region where the index value in the anisotropic region exceeds a threshold, and a maximum eigenvector having a maximum eigenvalue is determined from eigenvectors for each pixel. From the start pixel, the vector image of the maximum eigenvector is used to sequentially extract fibrous images from the start pixel.Therefore, the maximum eigenvector indicating the traveling direction of the nerve fiber is calculated, and the nerve fiber pixels along the direction of this vector are calculated. Can be sequentially and automatically extracted, and as a result, the time required for extracting nerve fiber pixels can be significantly reduced, and individual differences in extraction can be reduced. Camphor can be.
[Brief description of the drawings]
FIG. 1 is a block diagram illustrating an overall configuration of a magnetic resonance imaging system.
FIG. 2 is a diagram illustrating an imaging region and a pulse sequence according to the embodiment;
FIG. 3 is a diagram illustrating diffusion tensor image information according to the embodiment;
FIG. 4 is a functional block diagram of the image processing apparatus according to the embodiment;
FIG. 5 is a diagram illustrating diagonalization of a diffusion tensor matrix.
FIG. 6 is a diagram showing an operation of an area expanding unit.
FIG. 7 is a flowchart illustrating an operation of the image processing apparatus according to the embodiment;
FIG. 8 is a diagram showing the magnitude and direction of a diffusion coefficient in a nerve fiber.
FIG. 9 is a diagram schematically illustrating extraction of a nerve fiber image.
[Explanation of symbols]
1) Subject
20mm cabinet
30 operation console
40 image processing device
100mm magnet part
102 main magnetic field coil
106 ° gradient coil
108 RF coil section
120 cradle
130 ° gradient drive
140 ° RF drive
150 Data collection unit
155 Data processing unit
160 scan controller
170 ° control unit
171 Image processing unit
166, 176, 177 storage unit
180, 181 display
190, 191 operation section
200 diffusion tensor image information
202 T₂Enhanced image information
210 tensor matrix calculation means
220 matrix diagonalization means
230 ° index calculation means
240 ° threshold means
250 area expansion means
260 ° volume rendering means

Claims

Run the imaging pulse sequence to acquire the diffusion tensor image information and T ₂ weighted images information of the same imaging region of the object,
Calculating a diffusion tensor matrix for each pixel from the diffusion tensor image information and the T ₂ weighted images information,
Calculating eigenvalues and eigenvectors of the diffusion tensor matrix,
Calculating an index value of anisotropy of a diffusion coefficient at the pixel position from the eigenvalue,
Based on the index value, extract an anisotropic region from the imaging region,
A starting pixel is selected in a region where the index value in the anisotropic region exceeds a threshold,
From the eigenvectors for each pixel, select the largest eigenvector having the largest eigenvalue,
From the start pixel, using the vector information of the maximum eigenvector, to sequentially extract a fibrous image,
A fibrous image extraction method characterized by the above-mentioned.

Storage means for storing the diffusion tensor image information and T ₂ weighted images information of the same imaging region of the object,
And tensor matrix calculating means for calculating a diffusion tensor matrix for each pixel from the diffusion tensor image information and the T ₂ weighted images information,
Matrix diagonalization means for calculating eigenvalues and eigenvectors of the diffusion tensor matrix,
An index calculation unit that calculates an index value of anisotropy of a diffusion coefficient at the pixel position from the eigenvalue,
Based on the index value, a threshold means for extracting a high anisotropic region from the imaging site,
An image processing apparatus comprising:
A starting pixel is selected in a region where the index value in the anisotropic region exceeds a threshold,
From the eigenvectors for each pixel, select the largest eigenvector having the largest eigenvalue,
From the start pixel, using the vector information of the maximum eigenvector, region expansion means for extracting a fibrous image,
An image processing apparatus comprising:

The image processing apparatus according to claim 2, wherein the index value is a fractional anisotropy value calculated from the eigenvalue.

The image processing apparatus according to claim 3, wherein the threshold unit extracts a pixel region in which the fractional anisotropy value exceeds a threshold.

The image processing apparatus according to claim 2, wherein the tensor matrix calculation unit is optimized using a least squares method.

The area expansion unit calculates a connection index indicating a connection between a central pixel and a peripheral pixel of the mask area for each of the peripheral pixels from the vector information, and sets the peripheral pixel whose coupling index does not exceed a threshold value to a new center. The image processing apparatus according to any one of claims 2 to 5, wherein a new mask area as a pixel is set, and the area is sequentially expanded.

The image processing apparatus according to claim 6, wherein the connection index is an angle of a direction cosine calculated from an inner product of the vector information of the center pixel and the peripheral pixels.

Forming a static magnetic field, and a gradient magnetic field, transmits or receives RF signals, and controls the gradient magnetic fields and the RF signals, a diffusion tensor image information and T ₂ weighted images information of the same imaging region of the subject A magnetic resonance imaging apparatus to acquire;
Transfer means for transferring the diffusion tensor image information and T ₂ weighted images and information on the magnetic resonance imaging apparatus,
An image processing apparatus that performs image processing on the transferred diffusion tensor image information and T ₂ -weighted images information by said transfer means,
A magnetic resonance imaging system comprising:
The image processing apparatus, the diffusion tensor image information and the tensor matrix calculating means for calculating a diffusion tensor matrix for each pixel from the T ₂ -weighted images information, matrix diagonalization means for calculating the eigenvalues and eigenvectors of the diffusion tensor matrix Index calculating means for calculating an index value of anisotropy of a diffusion coefficient at the pixel position from the eigenvalue; threshold value means for extracting an anisotropic region from the imaging region based on the index value; The starting pixel is selected in a region where the index value in the anisotropic region exceeds a threshold, and among the eigenvectors for each of the pixels, the largest eigenvector having the largest eigenvalue is selected. Area expansion means for extracting a fibrous image using vector information,
A magnetic resonance imaging system comprising: