子宮の蠕動運動を可視化するための画像処理装置、画像処理方法および画像診断支援システム、ならびにそのためのコンピュータプログラムおよびそれを記録した記録媒体
【課題】子宮の蠕動運動を分かりやすく可視化して、子宮の画像診断を支援する。
【解決手段】この画像処理装置は、複数枚の子宮ＭＲ画像に対応した複数枚の子宮ＭＲ画像データを処理する。画像処理演算部３０は、複数枚の子宮ＭＲ画像データに対して時間軸方向のフーリエ変換を行い、子宮ＭＲ画像中の各画素について複数の周波数成分の位相データを生成するフーリエ変換処理部３３と、少なくとも１つの周波数成分の位相データを用いて、子宮ＭＲ画像各部における当該周波数成分の所定方向に関する位相勾配を表す位相勾配データを生成する位相勾配データ生成処理部４０と、この生成された位相勾配データに基づき、画像各部を位相勾配に応じて着色した位相勾配分布画像を表す位相勾配分布画像データを、子宮の蠕動を表す画像データ（蠕動の方向および／または速さを表すもの）として生成する位相勾配分布画像データ生成処理部４１とを備えている。
【解決手段】この画像処理装置は、複数枚の子宮ＭＲ画像に対応した複数枚の子宮ＭＲ画像データを処理する。画像処理演算部３０は、複数枚の子宮ＭＲ画像データに対して時間軸方向のフーリエ変換を行い、子宮ＭＲ画像中の各画素について複数の周波数成分の位相データを生成するフーリエ変換処理部３３と、少なくとも１つの周波数成分の位相データを用いて、子宮ＭＲ画像各部における当該周波数成分の所定方向に関する位相勾配を表す位相勾配データを生成する位相勾配データ生成処理部４０と、この生成された位相勾配データに基づき、画像各部を位相勾配に応じて着色した位相勾配分布画像を表す位相勾配分布画像データを、子宮の蠕動を表す画像データ（蠕動の方向および／または速さを表すもの）として生成する位相勾配分布画像データ生成処理部４１とを備えている。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、時系列に従って子宮を撮影して得た複数枚の画像データに基づいて子宮の蠕動運動を解析して可視化するための画像処理装置、画像処理方法および画像診断支援システム、ならびにそのためのコンピュータプログラムおよびそれを記録した記録媒体に関する。
【背景技術】
【０００２】
婦人科領域の画像診断では、超音波画像、Ｘ線撮影画像およびＭＲ（Magnetic Resonance; 磁気共鳴) 画像が用いられている。超音波画像には死角が存在し、分解能が低く、組織特異性が描画されないという欠点があり、また、Ｘ線撮影にはＸ線被曝の問題がある。これに対して、ＭＲ画像を用いた診断は、Ｘ線を使用しないため、妊娠中の患者に対しても安全に撮影が行え、かつ、優れた組織コントラストにより、機能的変化および形態的変化を顕著に描画することができる。
【０００３】
従来型のＭＲＩ（Magnetic Resonance Imaging）は、撮影時間が長く、子宮内膜が波打つような反復した収縮運動や子宮内膜の外側の筋層の反復運動のような子宮蠕動運動を捉えることは不可能であった。そのため、画像診断において、子宮の蠕動は全く考慮されず、子宮はあたかも静止した臓器であるかのように扱われてきた。
しかし、近年になって、ＨＡＳＴＥ（half-Fourier-acquisition single-shot turbo spin-echo）等の高速撮影法が開発されるに至り、ＭＲ画像の撮影に要する時間が短縮され、これに伴い、子宮蠕動運動を直接描画する動態撮像が可能となってきた（非特許文献１、２）。
【０００４】
女性生殖器は、年齢、月経周期およびホルモン等の影響を受ける。そのため、子宮は、他の臓器とは異なり、正常時でもその形状は一定ではない。すなわち、子宮は、排卵期、月経期、黄体期に、それぞれ異なる性質を持った運動を行っていると考えられている。そこで、子宮の蠕動運動を解析することにより、月経、妊能に関する各患者の症状に適した治療法の発見につながることが期待されている。
【非特許文献１】植田隆史他著「骨盤部（女性）撮像に必要な解剖・生理機能」、INNERVISION(18・9)、２００３年
【非特許文献２】富樫かおり著「ＭＲＩによる婦人科画像診断」、日本放射線技術学会雑誌第５９巻第８号、２００３年８月、ｐｐ８９５−９０３
【発明の開示】
【発明が解決しようとする課題】
【０００５】
しかし、現在のＭＲ画像による子宮蠕動運動の分析は、生のＭＲ画像の肉眼による目視観察に依存している。ところが、生のＭＲ画像中に現れる子宮の蠕動は、子宮内膜に現れる微少な波状突起の移動や、子宮筋層部分内における微妙な陰影（低輝度部）の移動として観察されるに過ぎない。より具体的に説明すると、骨盤部を撮影したＭＲ画像中では、周囲の臓器の活発な運動が現れており、このような周囲の臓器の運動の影響を受けて子宮の形状も絶えず変動している。このような状況の画像中では、子宮の蠕動は、ノイズに近い微弱な運動として観察されるに過ぎず、熟練した専門医師による極めて注意深い観察を必要とする。
【０００６】
そのため、ＭＲ画像に基づく子宮蠕動運動の評価は主観的にならざるをえず、適正な診断には、専門医師であっても、相当の熟練を要するという問題がある。
そこで、この発明の目的は、子宮の蠕動運動を分かりやすく可視化することができ、これにより、子宮の画像診断を支援することができる画像処理装置、画像処理方法および画像診断支援システム、ならびにそのためのコンピュータプログラムおよびそれを記録した媒体を提供することである。
【課題を解決するための手段】
【０００７】
上記の目的を達成するための請求項１記載の発明は、子宮の蠕動運動を解析して可視化するための画像処理装置であって、子宮の蠕動の周期よりも短い周期（たとえば、２〜３秒間隔）で時系列に従って撮影した複数枚の子宮ＭＲ（磁気共鳴）画像に対応した複数枚の子宮ＭＲ画像データ（個々の画素の輝度データで構成されるもの。とくに、いわゆる矢状断面画像）を記憶するための画像データ記憶手段と、この画像データ記憶手段に記憶された複数枚の子宮ＭＲ画像データに対して時間軸方向のフーリエ変換を行い、子宮ＭＲ画像中の各画素について複数の周波数成分（輝度の時間変化の周波数成分）の位相データを生成するフーリエ変換手段と、このフーリエ変換手段によって生成された少なくとも１つの周波数成分の位相データを用いて、子宮ＭＲ画像各部における当該周波数成分の所定方向に関する位相勾配（当該所定方向に隣接する画素との位相データの差分。微分値）を算出し、その位相勾配を表す位相勾配データを生成する位相勾配データ生成手段と、この位相勾配データ生成手段によって生成された位相勾配データに基づき、画像各部を位相勾配に応じて着色した位相勾配分布画像を表す位相勾配分布画像データを、子宮の蠕動を表す画像データ（蠕動の方向および／または速さを表すもの）として生成する位相勾配分布画像データ生成手段とを含むことを特徴とする子宮蠕動可視化のための画像処理装置である。
【０００８】
子宮の蠕動は、子宮内膜における波状の微少突起のほぼ一定速度での周期的な移動と、子宮筋層における局部的な状態変化部位のほぼ一定速度での周期的な移動とに大別される。一般に、排卵期には、子宮内膜および子宮筋層における蠕動方向は、精子の輸送のために、子宮の頚部から体部へと向かう方向である。月経期には、子宮内膜および子宮筋層における蠕動方向は、月経血の排出のために、子宮体部から子宮頚部へと向かう方向である。さらに、黄体期（分泌期）では、妊卵の着床や初期の妊卵の保持のため、内膜は蠕動を行っていないと考えられている。
【０００９】
この発明では、時系列に従う複数枚の子宮ＭＲ画像データに対して時間軸方向のフーリエ変換が施される。これにより、子宮ＭＲ画像を構成する個々の画素について、画素値（通常は輝度値）の時間変化が複数の周波数成分に分解され、個々の周波数成分の位相データが生成される。前述のように、子宮の蠕動は、ほぼ一定速度での子宮内膜および子宮筋層の周期的な運動であるので、この蠕動に対応した周波数成分に着目することによって、子宮の蠕動を解析することができる。
【００１０】
そこで、この発明では、フーリエ変換によって生成された少なくとも１つの周波数成分の位相データが用いられ、子宮ＭＲ画像各部における当該周波数成分の所定方向に関する位相勾配を表す位相勾配データが算出される。子宮の蠕動に対応した周波数成分の位相データに関して所定方向の勾配（位相勾配）を算出すると、この位相勾配は、子宮各部の運動速度（蠕動の速さおよび方向）を表す。このような位相勾配データに応じた色データを画素値とする位相勾配分布画像を形成すると、子宮各部の運動（方向および／または速さ）を一目で把握することができる。この発明では、このような位相勾配分布画像を表す位相勾配分布画像データが生成される。
【００１１】
したがって、子宮の蠕動に対応した周波数成分の位相データに対応した位相勾配分布画像データを用いて表示装置に位相勾配分布画像を表示することにより、子宮蠕動を明瞭に可視化した画像を提供することができる。これにより、子宮に関する画像診断を効果的に支援することができる。
位相勾配を求める方向は、操作者が指定することとしてもよいし、所定の一定方向に予め設定しておいてもよい。
【００１２】
前記画像処理装置は、前記フーリエ変換手段によって生成された少なくとも１つの周波数成分の位相データを用いて、子宮ＭＲ画像各部の当該周波数成分の位相分布を表す位相分布データ（位相画像データ）を生成する位相分布データ生成手段（位相画像データ生成手段）を含んでいてもよい。この場合、前記位相勾配データ生成手段は、前記位相分布データ生成手段によって生成された位相分布データ（各画素の位相データ）に基づいて、子宮ＭＲ画像の各部の位相勾配を算出するものであってもよい。
【００１３】
請求項２記載の発明は、前記画像データ記憶手段に記憶された複数枚の子宮ＭＲ画像データに対して、子宮画像部分の位置および形状の変動を補償するための幾何学変換を施し、この幾何学変換処理後の子宮ＭＲ画像データを生成する幾何学変換手段をさらに含み、前記フーリエ変換手段は、前記幾何学変換手段による幾何学変換処理後の子宮ＭＲ画像データに対してフーリエ変換を実行するものであることを特徴とする請求項１記載の画像処理装置である。
【００１４】
ＭＲＩによって撮影される画像中で、子宮の位置および形状は時間経過とともに変動する。そこで、この発明では、複数枚の子宮ＭＲ画像データに対して、子宮画像部分の位置および形状を合わせるための幾何学変換を施し、その後に、フーリエ変換を行うようにしている。これにより、より確実に、子宮の蠕動に対応した画像の変化を抽出することができる。
【００１５】
請求項３記載の発明は、前記フーリエ変換手段は、前記位相データのほかに、子宮ＭＲ画像の各画素についての複数の周波数成分の振幅データを生成するものであり、このフーリエ変換手段が生成する振幅データに基づいて、少なくとも１つの周波数成分を選択する周波数成分選択手段をさらに含み、前記位相勾配データ生成手段は、前記周波数成分選択手段によって選択された周波数成分に関して子宮ＭＲ画像各部の位相勾配データを生成するものであることを特徴とする請求項１または２記載の画像処理装置である。
【００１６】
この構成では、フーリエ変換によって、子宮ＭＲ画像の各画素について、複数の周波数成分に関する位相データおよび振幅データが生成される。子宮ＭＲ画像中に見られる運動は、子宮の蠕動ばかりではなく、周辺の臓器からの影響による移動および変形などのノイズ的な運動も含まれている。フーリエ変換によって分解される複数の周波数成分のなかには、このようなノイズ的な運動を反映したものもある。
【００１７】
子宮を含む画像領域に注目すると、子宮蠕動に対応した周波数成分は、ノイズ的な運動に対応した周波数成分に比較すると、大きな振幅を持つ。そこで、フーリエ変換によって得られる個々の周波数成分の振幅データに注目することによって、子宮蠕動を表す周波数成分を選択することができる。そして、このようにして選択された周波数成分に関して子宮ＭＲ画像の各部の位相勾配データを求め、その分布を表す位相勾配分布画像データを作成すれば、この位相勾配分布画像データに対応する疑似カラー画像を表示装置に表示して、子宮の蠕動を分かりやすく可視化することができる。
【００１８】
なお、前記位相分布データ生成手段（位相画像データ生成手段）が備えられる場合には、この位相分布データ生成手段は、前記周波数成分選択手段によって選択された周波数成分に関する位相分布データを生成するものであることが好ましい。
請求項４記載の発明は、前記周波数成分選択手段は、子宮ＭＲ画像の各画素について、最大の振幅データに対応した最大振幅周波数成分を抽出する最大振幅周波数抽出手段と、この最大振幅周波数抽出手段によって個々の画素について抽出された最大振幅周波数成分の子宮ＭＲ画像内での分布を表す最大振幅周波数分布画像データ（たとえば、最大周波数成分に対応した色データを画素値とする画像データ）を生成する最大振幅周波数分布画像データ生成手段とを含むことを特徴とする請求項３記載の画像処理装置である。
【００１９】
この構成によれば、最大振幅周波数分布画像データに基づいて、表示装置に最大振幅周波数分布画像（たとえば、最大周波数成分の分布を表す疑似カラー画像）を表示すれば、子宮ＭＲ画像中での最大振幅周波数成分の分布を視覚的に把握することができる。この最大振幅周波数分布画像に基づいて、子宮の蠕動に対応した周波数成分を特定することができる。
【００２０】
より具体的には、請求項５に記載されているように、前記周波数成分選択手段は、さらに、操作者による周波数成分選択入力を受け付ける周波数成分選択入力受け付け手段を含むものであってもよい。
この画像処理装置を適用した画像診断支援システムでは、たとえば、操作者による入力操作が可能な入力装置と、前記最大振幅周波数分布画像データ生成手段が生成する最大振幅周波数分布画像データを画像表示する表示装置とが備えられる。操作者が、表示装置に表示された最大振幅周波数分布画像から判断して、入力装置からいずれかの特定の周波数成分の選択入力を行うことにより、この入力が前記周波数成分選択入力受け付け手段によって受け付けられ、その周波数成分が選択されることになる。
【００２１】
さらに具体的に説明すると、最大振幅周波数分布画像中において、子宮ＭＲ画像中の子宮の一部または全部に近似した形状を形成する周波数成分が存在すれば、このような周波数成分は、子宮蠕動に対応した周波数成分である可能性が高い。そこで、操作者は、このような周波数成分を選択して入力装置から入力することになる。
医療の現場における画像診断では、完全に自動化されたシステムよりも、操作者（一般的には専門医）による手動操作の余地を残すことが好まれる場合が多い。これは、診断の基礎となる画像は、個々の患者に依存する多様な画像であり、また、患者の状態によっても異なるため、画像の自動解析が必ずしも成功するとは限らず、必要に応じて専門医の判断に委ねるステップを設ける方が、信頼性の高い画像診断が可能であり、かつ、診断結果に対する安心感を得やすいからである。さらにまた、画像解析の一部のステップを操作者の判断に委ねることで、画像処理装置の演算負荷が著しく軽減されるという効果も得られる。
【００２２】
ただし、子宮蠕動に対応した周波数成分の選択の自動化が不可能であるというわけではなく、必要に応じて、周波数成分選択処理の自動化も可能である。
具体的には、請求項６に記載されているように、前記周波数成分選択手段は、さらに、前記最大振幅周波数分布画像データ生成手段が生成する最大振幅周波数分布画像データに基づき、当該最大振幅周波数分布画像データに対応した画像中に占める個々の周波数成分の面積の順位に基づいて、少なくとも１つの周波数成分を選択する面積順位基準選択手段を含むものであってもよい。
【００２３】
最大振幅周波数分布画像中では、子宮蠕動に対応した周波数成分の領域が支配的になる可能性が高い。したがって、最大振幅周波数分布画像中に占める面積が上位（たとえば第１位のみ。第１位〜第３位など）の周波数成分は、子宮蠕動に対応した周波数成分である可能性が高い。したがって、最大振幅周波数分布画像中での占有面積が上位の周波数成分を自動選択する前記の構成によって、操作者による判断を要することなく、子宮蠕動に対応した周波数成分を適切に選択することができる。
【００２４】
この構成は、フーリエ変換以降の処理対象となる子宮ＭＲ画像が、子宮近傍の領域を切り出した画像である場合にとくに効果的である。
また、請求項７に記載されているように、前記周波数成分選択手段が、所定のしきい値以上の最大振幅データに対応した周波数成分を選択するものである場合も、子宮蠕動に対応した周波数成分の自動選択が可能である。すなわち、子宮蠕動に対応した周波数成分はある程度の振幅を示すと考えられるので、所定のしきい値未満の最大振幅データに対応した周波数成分はノイズとみなせる。また、一定以上の振幅を持たない周波数成分については、その位相データに基づく解析自体があまり意味を持たなくなるから、無意味な画像解析処理を省くためにも、最大振幅データが前記しきい値に達しない周波数成分は、処理の対象外とすることが好ましい。
【００２５】
請求項８記載の発明は、前記フーリエ変換手段は、前記位相データのほかに、子宮ＭＲ画像の各画素についての複数の周波数成分の振幅データを生成するものであり、前記フーリエ変換手段によって生成された少なくとも１つの周波数成分の振幅データを用いて、子宮ＭＲ画像各部の当該周波数成分の振幅分布を表す振幅画像を表す振幅画像データを生成する振幅画像データ生成手段をさらに含むことを特徴とする請求項１ないし７のいずれかに記載の画像処理装置である。
【００２６】
この構成によれば、振幅画像データ生成手段によって生成される振幅画像データに基づいて表示装置における画像表示を行うことで、操作者に対して、個々の周波数成分の振幅分布を表す振幅画像を提供できる。この振幅画像は、たとえば、当該周波数成分の振幅データに対応する輝度を各画素に与えたものであってもよい。子宮蠕動に対応した周波数成分の振幅画像では、振幅の大きな画素が子宮の一部または全部の形状を形成する可能性が高い。したがって、このような振幅画像を提供することによって、操作者は、子宮蠕動に対応した周波数成分を容易に特定して選択することができる。
【００２７】
請求項９記載の発明は、前記フーリエ変換手段によって生成された少なくとも１つの周波数成分の位相データを用いて、子宮ＭＲ画像各部の当該周波数成分の位相分布を表す位相画像を表す位相画像データを生成する位相画像データ生成手段をさらに含むことを特徴とする請求項１ないし８のいずれかに記載の画像処理装置である。
この構成によれば、位相画像データ生成手段によって生成される位相画像データに基づいて表示装置における画像表示を行うことで、操作者に対して、個々の周波数成分の位相分布を表す位相画像を提供することができる。たとえば、この位相画像は、子宮ＭＲ画像の当該周波数成分の位相に対応した輝度を個々の画素に与えたものであってもよい。
【００２８】
子宮蠕動に対応した周波数成分について前述のような位相画像を作成すると、この位相画像には、子宮の一部または全部の形状に従って位相変化を示す画像部分が現れる。したがって、位相画像を提供することによって、子宮蠕動に対応する周波数成分を特定して選択する過程で、操作者を支援することができる。
或る周波数成分の位相画像を作成するときには、当該周波数成分の振幅が所定のしきい値に満たない画素については、画素値を固定（たとえば「０」に固定）することが好ましい。これにより、子宮蠕動に全く関係のないと考えられる部分の画像を潰すことができるので、操作者による判断がより容易になる。
【００２９】
請求項１０記載の発明は、前記画像データ記憶手段に記憶された子宮ＭＲ画像データ（複数枚の子宮ＭＲ画像データのなかから選択した代表的な１枚であってもよいし、複数の子宮ＭＲ画像データの平均値等をとった画像データであってもよい。）に前記位相勾配分布画像データ生成手段によって生成される位相勾配分布画像データをオーバレイして、オーバレイ画像データを生成する画像オーバレイ処理手段をさらに含むことを特徴とする請求項１ないし９のいずれかに記載の画像処理装置である。
【００３０】
この構成によれば、オーバレイ画像データに基づいて表示装置上に画像表示を行えば、子宮ＭＲ画像に位相勾配分布を表す着色画像がオーバレイ表示される。これにより、子宮の各部と、蠕動が生じている部位との対応関係を一目で把握できる画像を提供することができ、どの部位でどのような蠕動が生じているかを容易に判断できるようになる。
請求項１１記載の発明は、前記フーリエ変換手段によって生成された１つの周波数成分の位相データのうち、少なくとも２つの離間した画素の位相データに基づいて、当該２つの画素間の位相差を演算する位相差演算手段と、この位相差演算手段によって演算された位相差に基づいて、子宮の蠕動速度を計算する蠕動速度演算手段をさらに含むことを特徴とする請求項１ないし１０のいずれかに記載の画像処理装置である。
【００３１】
子宮蠕動に対応した周波数成分に関して２つの画素間の位相差を求めると、この位相差に基づいて子宮蠕動の速度を算出することができる。この算出された子宮蠕動の速度を、たとえば、表示装置に表示したりして提供することにより、子宮蠕動に関するより多くの情報を提供でき、より一層効果的に画像診断を支援できる。子宮蠕動の速度は、より具体的には、２つの画素間の位相差と、この２つの画素間の距離とに基づいて、算出することができる。２つの画素間の距離は、２つの画素の位置情報（座標）から算出できる。
【００３２】
請求項１２記載の発明は、操作者による画像中の２点の指示入力を受け付ける手段をさらに含み、前記位相差演算手段は、前記受け付けられた指示入力による２点の画素の位相データ（同じ周波数成分の位相データ）に基づいて、前記２点の画素間の位相差を演算するものであることを特徴とする請求項１１記載の画像処理装置である。
この構成によれば、入力装置から、速度の計測を行いたい２点を指示することができるので、操作者の所望する任意の２点間の蠕動速度の情報を提供することができる。
【００３３】
請求項１３記載の発明は、請求項１ないし１２のいずれかに記載の画像処理装置と、この画像処理装置に接続され、前記位相勾配分布画像を表示する表示装置とを含むことを特徴とする画像診断支援システムである。
この構成により、位相勾配分布画像データに対応した位相勾配分布画像を表示装置に表示して、この画像情報を使用者に提供することができる。これにより、子宮の蠕動の様子を明瞭に可視化した画像を提供でき、画像診断を良好に支援できる。
【００３４】
前記表示装置には、前述の最大振幅周波数分布画像、振幅画像および位相画像の表示も可能とされていることが好ましい。
請求項１４記載の発明は、操作者によって操作可能であり、前記画像処理装置に接続されて、この画像処理装置に対して操作者の指示入力を入力するための入力装置をさらに含むことを特徴とする請求項１３記載の画像診断支援システムである。
【００３５】
この構成により、使用者による指示入力によって、画像処理装置による処理内容を指示することができる。
請求項１５記載の発明は、子宮の蠕動運動を解析して可視化するための画像処理方法であって、子宮の蠕動の周期よりも短い周期で時系列に従って撮影した複数枚の子宮ＭＲ（磁気共鳴）画像に対応した複数枚の子宮ＭＲ画像データを準備するステップと、この複数枚の子宮ＭＲ画像データに対して時間軸方向のフーリエ変換を行い、子宮ＭＲ画像中の各画素について複数の周波数成分の位相データを生成するフーリエ変換ステップと、このフーリエ変換ステップによって生成された少なくとも１つの周波数成分の位相データを用いて、子宮ＭＲ画像各部における当該周波数成分の所定方向に関する位相勾配を算出し、その位相勾配を表す位相勾配データを生成する位相勾配データ生成ステップと、この位相勾配データ生成ステップによって生成された位相勾配データに基づき、画像各部を位相勾配に応じて着色した位相勾配分布画像を表す位相勾配分布画像データを、子宮の蠕動を表す画像データとして生成する位相勾配分布画像データ生成ステップとを含むことを特徴とする子宮蠕動可視化のための画像処理方法である。この方法により、請求項１の発明に関して説明した効果と同様の効果が得られる。
【００３６】
請求項１６記載の発明は、前記フーリエ変換ステップによって生成された１つの周波数成分の位相データのうち、少なくとも２つの離間した画素の位相データに基づいて、当該２つの画素間の位相差を演算する位相差演算ステップと、この位相差演算ステップによって演算された位相差に基づいて、子宮の蠕動速度を計算する蠕動速度演算ステップとをさらに含むことを特徴とする請求項１５記載の画像処理方法である。この方法により、請求項１１の発明と同様な効果が得られる。
【００３７】
むろん、画像処理方法の発明に関しても、画像処理装置の発明と同様な変形および改良が可能である。
請求項１７記載の発明は、子宮の蠕動の周期よりも短い周期で時系列に従って撮影した複数枚の子宮ＭＲ（磁気共鳴）画像に対応した複数枚の子宮ＭＲ画像データを処理し、子宮の蠕動運動を解析して可視化するための画像処理装置としてコンピュータを動作させるためのコンピュータプログラムであって、前記コンピュータを、前記複数枚の子宮ＭＲ画像データに対して時間軸方向のフーリエ変換を行い、子宮ＭＲ画像中の各画素について複数の周波数成分の位相データを生成するフーリエ変換手段、このフーリエ変換手段によって生成された少なくとも１つの周波数成分の位相データを用いて、子宮ＭＲ画像各部における当該周波数成分の所定方向に関する位相勾配を算出し、その位相勾配を表す位相勾配データを生成する位相勾配データ生成手段、およびこの位相勾配データ生成手段によって生成された位相勾配データに基づき、画像各部を位相勾配に応じて着色した位相勾配分布画像を表す位相勾配分布画像データを、子宮の蠕動を表す画像データとして生成する位相勾配分布画像データ生成手段として機能させることを特徴とする子宮蠕動可視化のためのコンピュータプログラムである。このコンピュータプログラムをコンピュータによって実行させることにより、このコンピュータを請求項１記載の画像処理装置として機能させることができる。
【００３８】
請求項１８記載の発明は、前記コンピュータを、さらに、前記フーリエ変換手段によって生成された１つの周波数成分の位相データのうち、少なくとも２つの離間した画素の位相データに基づいて、当該２つの画素間の位相差を演算する位相差演算手段、およびこの位相差演算手段によって演算された位相差に基づいて、子宮の蠕動速度を計算する蠕動速度演算手段として機能させることを特徴とする請求項１７記載のコンピュータプログラムである。このコンピュータプログラムをコンピュータに実行させることにより、このコンピュータを請求項１１の画像処理装置として機能させることができる。
【００３９】
むろん、コンピュータプログラムは、画像処理装置の発明に関連して説明した前述の種々の機能処理手段としてコンピュータを機能させるように変形または改良されてもよい。
請求項１９記載の発明は、前述のようなコンピュータプログラムを記録しており、コンピュータによる読み取りが可能な記録媒体である。このような記録媒体としては、ＣＤ−ＲＯＭおよびＤＶＤ−ＲＯＭに代表される光ディスク、ＭＯディスクに代表される磁気光ディスク、ならびにハードディスクドライブおよびフレキシブルディスクに代表される磁気ディスクを例示することができる。
【発明を実施するための最良の形態】
【００４０】
以下では、この発明の実施の形態を、添付図面を参照して詳細に説明する。
［子宮蠕動の解析の基本原理］
まず、この発明の一実施形態に係る画像処理装置による子宮ＭＲ画像の解析の基本原理について説明する。
図１に、子宮ＭＲ画像（矢状断面）の一例を示す。子宮は、領域ａに示す「体部」と、領域ｂに示す頚部とを有している。子宮の壁は、３層で構成されている。一番内側が内膜ｃであり、その内膜を覆うように存在するのが筋層ｄであり、内膜ｃと筋層ｄの境界がジャンクショナルゾーン（junctional zone）である。
【００４１】
図１に示す子宮ＭＲ画像は、高速スピンエコー法によって撮影された画像例である。スピンエコー法では、Ｔ１強調画像およびＴ２強調画像の撮影が可能であるが、Ｔ１強調画像では、子宮筋層と内膜が均一な低信号を持つため、これらを区別して撮影することができない。そのため、子宮のＭＲＩでは、Ｔ２強調画像での撮影が主として行われ、図１に示す子宮ＭＲ画像もＴ２強調画像である。
【００４２】
子宮の蠕動の解析には、子宮の蠕動の周期よりも短い時間間隔（たとえば２秒間隔）で時系列に従って撮影された複数枚の子宮ＭＲ画像のデータが用いられる。画像サイズは、たとえば、横１．２ｍｍ、縦１．４ｍｍである。撮影の性質上、１枚目の画像の輝度値は、他の画像に比べて高くなるので、１枚目は除き、２枚目からの画像を使用して解析を行うことが好ましい。後述する高速フーリエ変換のためには、２のべき乗枚（たとえば６４枚）の子宮ＭＲ画像が必要となる。したがって、たとえば、６０枚の子宮ＭＲ画像データが処理対象であり、２枚目以降の５９枚の子宮ＭＲ画像データに対して画像解析を行う場合には、これを、補間処理によって、６４枚の子宮ＭＲ画像データに変換して使用することになる。
【００４３】
複数枚の子宮ＭＲ画像をディスプレイ上に連続で表示して注意深く観察すると、子宮の内膜に内側に沿って微小な波状突起がほぼ一定の速さで移動し、子宮の内膜が波打つように運動をしていることが観察される。この微小波状突起の移動に伴い、内膜の画像部分を構成する各画素の輝度値が周期的な時間変化を示す。また、筋層では、黒い影のようなものが動いていることが観察できる。筋層は厚い筋肉から成り立っており、内膜と同様、周期的な運動をしていると考えられる。
【００４４】
図２に示す子宮ＭＲ画像において、子宮内膜上の点「１」、「２」および「３」の３点について、時間軸方向の輝度値の変化をプロットした例を図３に示す。横軸は、時系列に従う複数枚の画像の番号によって表される時間軸である。この図３から、子宮内膜の蠕動運動は、定点に注目したときに、濃度（輝度）変化の周期的な波の組合せとして表現できると推測できる。子宮筋層の蠕動運動についても同様である。
【００４５】
そこで、子宮ＭＲ画像の各画素の時系列データ（時系列に従って撮影された複数枚の子宮ＭＲ画像における同一位置の画素の輝度データ列）に対して、フーリエ変換処理が行われる。実際には、図２に示すように子宮部分以外の画像部分を排除するように矩形領域にトリミングした子宮ＭＲ画像を対象画像としてフーリエ変換処理が行われる。
ただし、子宮の画像部分の位置やその形状は、周囲の臓器の運動等の影響を受けて変動する。そのため、処理領域を固定して解析を行うフーリエ変換では、子宮の位置が変化すると、運動解析の対象位置がずれてしまい、適切な解析を行うことができない。そこで、フーリエ変換に先だって、複数枚の子宮ＭＲ画像間の子宮の位置および形状の変動を補償するための幾何学変換処理が行われる。具体的には、画像間の変形を近似するアフィン変換パラメータを非線形最適化を用いて求め、そのパラメータを用いたアフィン変換によって、子宮の位置および形状が合わせ込まれる。より詳細には、まず並進移動成分を決定した後、その並進移動のパラメータを初期値としてアフィン変換パラメータが求められる。このアフィン変換パラメータを用いて原画像が変形され、対応する領域の切り出し（トリミング）を行った画像列を作成し、その画像列に対してフーリエ変換が行われる。
【００４６】
フーリエ変換には、離散フーリエ変換の一形態である高速フーリエ変換（ＦＦＴ：Fast Fourier Transform）が適用される。これにより、各画素の輝度変化における主な周波数成分の振幅および位相が求められる。これに基づいて、子宮内膜および筋層の蠕動の方向および速度を解析することができる。
より具体的に説明すると、Ｎ枚の子宮ＭＲ画像中の同じ位置にある１つの画素についてのＮ個の時系列データ（輝度値を表す画像データ）をｘ0，ｘ1，ｘ2，…，ｘn，…，ｘN-1と表すと、周波数ωk＝２πｋ／Ｎ（ｋ＝０，１，２，３，……，Ｎ−１）についての離散フーリエ変換Ｘkは、次の式(1)で与えられる。ただし、ｊは虚数単位である。
【００４７】
【数１】
【００４８】
そして、周波数ωkの周波数成分の振幅データ｜Ｃk｜および位相データφkは、それぞれ、次の式(2)および式(3)で与えられる。
｜Ｃk｜＝｛ＸkRe2＋ＸkIm2｝1/2 ……（２）
φk ＝ｔａｎ-1（−ＸkIm／ＸkRe） ……（３）
ただし、ＸkReおよびＸkImは、それぞれ、Ｘkの実数部および虚数部である。
【００４９】
これらの振幅データ｜Ｃk｜および位相データφkが、子宮ＭＲ画像の個々の画素について、計算可能なすべての周波数ωkの周波数成分に関して求められる。
高速フーリエ変換を行うためには、前述のように、２のべき乗個の時系列データが必要となる。したがって、有効な子宮ＭＲ画像が５９枚であるときは、この５９枚の子宮ＭＲ画像（アフィン変換後に切り出した画像）の個々の画素の輝度値を５９個から６４個に補間したデータを予め作成し、この６４個のデータに対して高速フーリエ変換が行われる。
【００５０】
この高速フーリエ変換によって得られた振幅データに基づき、周波数「０」（直流成分）以外で最大振幅をとる周波数が、各画素について求められる。そして、各画素の最大振幅の周波数成分を表す最大振幅周波数分布画像がディプレイ上に表示される。
具体的には、図４(a)に示すように、複数の周波数成分に対して予め異なる色が割り当てられ、個々の画素が、当該画素の最大振幅の周波数成分に応じて着色される。こうして、最大振幅の周波数成分分布が、ディプレイ上に疑似カラー表示されることになる。
【００５１】
子宮の蠕動運動は低周波成分に現れるため、図４(a)の例では、第１周波数〜第１０周波数（ｋ＝１〜１０）の周波数成分についての分布画像を求めた例が示されている。第１周波数とは、６４枚で１周期となる周波数成分であり、第２周波数とは、３２（＝６４／２）枚で１周期となる周波数成分である。以下、第３周波数成分以上の周波数成分についても同様である。
【００５２】
参考のため、図４(b)には図４(a)のカラー画像中の第１周波数の色部分を強調した画像を示し、図４(c)には図４(a)のカラー画像中の第２周波数の色部分を強調した画像を示し、図４(d)には図４(a)のカラー画像中の第６周波数の色部分を強調した画像を示す。
内膜部分または筋層部分の蠕動は、規則的な周期運動であるので、内膜部分または筋層部分の画素値（最大振幅の周波数成分）は同じであることが期待される。実際、図４の例では、第６周波数成分を最大振幅周波数成分とする複数の画素が、子宮内膜に近似した形状の画像部分を形成している。
【００５３】
したがって、最大振幅周波数分布画像を見ることによって、使用者（操作者。一般的には専門医）は、子宮の蠕動に対応した周波数成分を容易に特定することができる。
子宮蠕動に対応した周波数成分の特定を補助する目的で、必要に応じて、振幅画像データおよび位相画像データが作成され、それらに対応する振幅画像および位相画像がそれぞれディスプレイ上に表示される。
【００５４】
振幅画像データは、一つの周波数成分について、各画素の振幅データを画素値（輝度値）とした振幅分布データである。このような振幅画像データに対応する振幅画像は、当該周波数成分の振幅が大きい画素ほど明るく表示される。子宮の蠕動は、前述のとおり、規則的な周期運動であるから、子宮の蠕動に対応した周波数成分についての振幅画像は、子宮の内膜部分または筋層部分に近似した形状の画像部分を形成すると期待される。したがって、振幅画像中に子宮の形状に近似した画像部分が形成される周波数成分を、子宮蠕動に対応するものと特定することができる。振幅画像の一例は、図１６に示されている。
【００５５】
位相画像データは、一つの周波数成分について、各画素の位相データを画素値（輝度）とした位相分布データである。ただし、当該周波数成分についての振幅データが所定のしきい値未満の画素については、画素値を「０」に固定する。このような位相画像データは、位相が大きい画素ほど明るく表示される。しきい値未満の画素の画素値を「０」とするので、子宮蠕動に対応した周波数成分の位相画像中には、やはり、子宮の形状に対応した画像部分が形成されることが期待でき、さらに、この画像部分は、子宮の蠕動方向および速度に対応した輝度変化を示すと期待される。したがって、このような特徴の画像部分が位相画像中に見られる周波数成分を、子宮蠕動に対応するものと特定することができる。位相画像の一例は、図１７に示されている。
【００５６】
こうして子宮蠕動に対応した周波数成分が特定されると、位相画像の濃度勾配（輝度勾配）が算出される。すなわち、個々の画素について、所定方向に関して、隣接画素との画素値の差分（濃度勾配）が求められる。この差分は、結局、位相勾配を表す位相勾配データにほかならない。位相勾配を求める方向は、たとえば、操作者がディスプレイに表示された位相画像上で指定する。
【００５７】
さらに、位相勾配データの値に応じて色情報が割り当てられる。たとえば、位相勾配データが正であれば、第１の色（たとえば青色）が当該画素に割り当てられ、位相勾配データが負であれば、第１の色とは異なる第２の色（たとえば赤色）が当該画素に割り当てられる。これにより、子宮蠕動の方向を色情報で表した位相勾配分布画像が形成される。位相勾配分布画像の例は、図２２、図２８および図３６に示されている。
【００５８】
むろん、位相勾配の値に応じて画素値（輝度）を異ならせてもよい。これにより、子宮蠕動の速さを併せて表示することができる。
子宮蠕動の速度は、画像中の任意の２点を指定することによって行える。速度＝距離÷時間であり、時間＝周期×位相差÷２πである。「周期」は１周期の時間であり、周期＝画像枚数×２（秒）÷周波数で表される。周波数とは、子宮蠕動に対応すると特定された周波数成分の周波数である。
【００５９】
したがって、子宮蠕動に対応すると特定された周波数成分について、指定された２点の位相データを用いて位相差を演算すれば、この位相差に基づいて、子宮蠕動の速度を求めることができる。
［画像診断支援システムの構成］
図５は、この発明の一実施形態に係る画像診断支援システムの基本的なハードウェア構成を説明するためのブロック図である。この画像診断支援システムは、コンピュータシステムとしての基本構成を有している。すなわち、この画像診断支援システムは、コンピュータ本体１０と、このコンピュータ本体１０に接続された入力装置１１としてのキーボード１１Ａおよびマウス（ポインティングデバイスの一例）１１Ｂと、コンピュータ本体１０に接続されたディプレイ（表示装置）１３と、コンピュータ本体１０に接続されたプリンタ１４とを備えている。
【００６０】
コンピュータ本体１０は、演算処理部としてのＣＰＵ１５と、主記憶装置としてのＲＯＭ１６およびＲＡＭ１７と、補助記憶装置してのハードディスクドライブ（ＨＤＤ）１８およびＣＤ−ＲＯＭドライブ１９と、入力制御部２０と、ディスプレイインタフェース２１と、プリンタインタフェース２２とを備えており、これらはバス２３によって互いに接続されている。入力制御部２０は、入力装置１１からの入力信号を処理してＣＰＵ１５に入力する。ディスプレイインタフェース２１は、画像データの入力を受けて、その画像データに対応する画像をディスプレイ１３に表示させる働きを担う。プリンタインタフェース２２は、プリンタ１４との接続のためのものである。
【００６１】
ＣＤ−ＲＯＭドライブ１９には、コンピュータプログラムを記録した記録媒体としてのＣＤ−ＲＯＭ（ディスク）２５を装填することができる。画像処理プログラムを記録したＣＤ−ＲＯＭ２５をＣＤ−ＲＯＭドライブ１９に装填し、必要なインストール操作を行ったうえで、当該画像処理プログラムをコンピュータ本体１０によって実行させることにより、このコンピュータ本体１０は画像処理装置として機能し、当該コンピュータシステムは画像診断支援システムとして機能することができる。
【００６２】
むろん、コンピュータプログラムは、記録媒体に記録された形態で提供される必要はなく、ネットワークを介して提供されてもよい。すなわち、コンピュータ本体１０にネットワークインタフェースを備え、このネットワークインタフェースを介してコンピュータプログラムをコンピュータ本体１０にダウンロードすることによっても、コンピュータ本体１０を画像処理装置として機能させることができる。
【００６３】
一方、処理対象の子宮ＭＲ画像データは、たとえば、ハードディスクドライブ１８に格納される。この子宮ＭＲ画像データも、ＣＤ−ＲＯＭドライブ１９を経由してハードディスクドライブ１８に格納されてもよいし、ネットワークを介して取得してハードディスクドライブ１８に格納するようにしてもよい。子宮ＭＲ画像データは、子宮の蠕動運動の周期よりも短い時間間隔（たとえば２秒）で時系列に従って撮影された複数枚の子宮ＭＲ画像を表すデータである。
【００６４】
図６は、前記画像処理装置として機能するコンピュータ本体１０の機能的な構成を説明するためのブロック図である。コンピュータ本体１０は、ＣＰＵ１５によって画像処理プログラムを実行させることによって、図６に示す機能処理部を実質的に有することになる。すなわち、コンピュータ本体１０は、画像処理演算部３０と、データ記憶部５０と、画像表示制御部７０とを有している。画像処理演算部３０および画像表示制御部７０は、主として、ＣＰＵ１５による演算処理によって実現される機能処理部である。また、データ記憶部５０は、ＲＡＭ１７およびハードディスクドライブ１８の記憶領域で構成される。
【００６５】
画像処理演算部３０は、アフィン変換処理部３１と、トリミング処理部３２と、フーリエ変換処理部３３と、最大振幅周波数成分抽出処理部３４と、最大振幅周波数分布画像データ生成処理部３５と、周波数成分選択入力受付け部３６と、振幅画像データ生成処理部３７と、位相画像データ生成処理部３８と、位相勾配算出方向入力受付け部３９と、位相勾配データ生成処理部４０と、位相勾配分布画像データ生成処理部４１と、画像オーバレイ処理部４２と、速度測定点入力受付け部４３と、位相差演算部４４と、蠕動速度演算部４５とを備えている。
【００６６】
また、データ記憶部５０は、原子宮ＭＲ画像データ記憶部５１と、アフィン変換後子宮ＭＲ画像データ記憶部５２と、対象子宮ＭＲ画像データ記憶部５３と、フーリエ変換データ記憶部５４と、最大振幅周波数成分データ記憶部５５と、最大振幅周波数分布画像データ記憶部５６と、振幅画像データ記憶部５７と、位相画像データ記憶部５８と、位相勾配データ記憶部５９と、位相勾配分布画像データ記憶部６０と、オーバレイ画像データ記憶部６１と、位相差データ記憶部６２と、蠕動速度データ記憶部６３とを備えている。
【００６７】
原子宮ＭＲ画像データ記憶部５１は、ＭＲＩによって子宮の蠕動運動の周期よりも短い時間間隔（たとえば２秒）で時系列に従って撮影された複数枚の子宮ＭＲ画像（原データ。ただし、必要に応じて補間処理をしたもの）を記憶する。
アフィン変換処理部３１は、原子宮ＭＲ画像データ記憶部５１から子宮ＭＲ画像データを読み出し、この子宮ＭＲ画像データにアフィン変換処理を行い、アフィン変換後の子宮ＭＲ画像データをアフィン変換後子宮ＭＲ画像データ記憶部５２に格納する。
【００６８】
トリミング処理部３２は、アフィン変換処理後の子宮ＭＲ画像データをアフィン変換後子宮ＭＲ画像データ記憶部５２から読み出し、この子宮ＭＲ画像データから子宮部分以外の画像を排除するトリミング処理を行い、この処理後のトリミングされた子宮ＭＲ画像データを対象子宮ＭＲ画像データ記憶部５３に格納する。
フーリエ変換処理部３３は、このトリミング処理部３２による処理後の子宮ＭＲ画像データを対象子宮ＭＲ画像データ記憶部５３から読み出し、この子宮ＭＲ画像データに対してフーリエ変換（高速フーリエ変換）を行う。より具体的には、フーリエ変換処理部３３は、アフィン変換を経てトリミングされた複数枚の子宮ＭＲ画像データに対して、時間軸方向のフーリエ変換を行う。これにより、子宮ＭＲ画像の個々の画素について、複数の周波数成分に関して、振幅データおよび位相データが求められる。すなわち、子宮ＭＲ画像の各画素の輝度値の時間変化が複数の周波数成分に分解され、各周波数成分について、その振幅を表す振幅データおよびその位相を表す位相データが求められる。これらのデータは、フーリエ変換データ記憶部５４に格納される。
【００６９】
最大振幅周波数成分抽出処理部３４は、フーリエ変換データ記憶部５４に格納された各画素の複数の周波数成分についての振幅データに基づき、子宮ＭＲ画像の個々の画素について、最大の振幅データを有する周波数成分を抽出する。こうして、画素ごとに抽出された周波数成分を表すデータは、最大振幅周波数成分データ記憶部５５に記憶される。
最大振幅周波数分布画像データ生成処理部３５は、最大振幅周波数成分データ記憶部５５に記憶されたデータに基づき、前記抽出された最大振幅周波数成分の分布を表す画像（最大振幅周波数分布画像）を表す最大振幅周波数分布画像データを生成し、最大振幅周波数分布画像データ記憶部５６に格納する。具体的には、最大振幅周波数分布画像データ生成処理部３５は、複数の周波数成分に対して異なる色データを割り当て、当該色データを画素値とする最大振幅周波数分布画像データを生成する。画像表示制御部７０は、この最大振幅周波数分布画像データを読み出して、ディスプレイ１３（図５参照）に最大振幅周波数分布画像を表示させる。これにより、子宮ＭＲ画像中における最大振幅周波数成分の分布を表す疑似カラー表示画像がディスプレイ１３上に表示される。
【００７０】
周波数成分選択入力受付け部３６は、入力装置１１からの周波数成分選択入力を受け付ける。操作者は、ディスプレイ１３に表示された最大振幅周波数分布画像から判断して、子宮蠕動に対応した周波数成分を特定し、この周波数成分の選択指示を入力装置１１を介して入力する。この選択指示が、周波数成分選択入力受付け部３６によって受け付けられる。
【００７１】
振幅画像データ生成処理部３７は、周波数成分選択入力受付け部３６によって受け付けられた選択指示に対応する周波数成分についての振幅画像の表示のための振幅画像データを生成する。具体的には、振幅画像データ生成処理部３７は、フーリエ変換データ記憶部５４から、前記選択された周波数成分についての各画素の振幅データを読み出し、この振幅データを画素値とする振幅画像データ（振幅分布データ）を作成し、この振幅画像データを振幅画像データ記憶部５７に格納する。画像表示制御部７０は、この振幅画像データを読み出して、ディスプレイ１３の表示制御を行う。これにより、前記選択された周波数成分についての振幅データを各画素の輝度値とする振幅画像がディスプレイ１３上に表示される。
【００７２】
一方、位相画像データ生成処理部３８は、周波数成分選択入力受付け部３６によって受け付けられた選択指示に対応する周波数成分についての位相画像の表示のための位相画像データを生成する。具体的には、位相画像データ生成処理部３８は、フーリエ変換データ記憶部５４から、前記選択された周波数成分についての各画素の位相データを読み出し、この位相データを画素値とする位相画像データ（位相分布データ）を作成し、この位相画像データを位相画像データ記憶部５８に格納する。ただし、位相画像データ生成処理部３８は、振幅画像データ記憶部５７に記憶された振幅画像データを参照し、その値（すなわち振幅データの値）が、所定のしきい値に達していない画素については、位相画像データを「０」に固定する。画像表示制御部７０は、位相画像データ記憶部５８から位相画像データを読み出して、ディスプレイ１３の表示制御を行う。これにより、位相画像データが「０」に固定された画素以外の画素について、前記選択された周波数成分についての位相データを各画素の輝度値とする位相画像がディスプレイ１３上に表示されることになる。
【００７３】
位相勾配算出方向入力受付け部３９は、入力装置１１からの位相勾配算出方向指示入力を受け付ける。操作者は、ディスプレイ１３に表示される振幅画像および／または位相画像から、子宮蠕動の方向を判断し、その方向を入力装置１１から指示する。この指示入力が位相勾配算出方向入力受付け部３９によって受け付けられることになる。
位相勾配データ生成処理部４０は、位相画像データ記憶部５８に記憶された位相画像データを読み出し、個々の画素について、前記位相勾配算出方向入力受付け部３９によって受け付けられた位相勾配算出方向に隣接する画素との間での差分値（微分値）を求め、この差分値を位相勾配データとして生成し、位相勾配データ記憶部５９に格納する。位相勾配データは、したがって、子宮ＭＲ画像の各部における運動の速度（方向および速さ）を表す。位相画像データ記憶部５８に、子宮蠕動に対応した周波数成分の位相画像データが格納されていると、位相勾配データは、子宮の各部の蠕動の速度（方向および速さ）を表すことになる。
【００７４】
位相勾配分布画像データ生成処理部４１は、位相勾配データ記憶部５９から位相勾配データを読み出し、位相勾配分布をカラー表示するための位相勾配分布画像データを生成する。具体的には、位相勾配分布画像データ生成処理部４１は、位相勾配データの値ごとに予め異なる色データを割り当てており、このような色データを画素値とする画像データを位相勾配分布画像データとして生成し、位相勾配分布画像データ記憶部６０に格納する。画像表示制御部７０は、位相勾配分布画像データ記憶部６０から位相勾配分布画像データを読み出し、ディスプレイ１３の表示制御処理を実行する。これにより、位相勾配分布を疑似カラー表示した位相勾配分布画像がディスプレイ１３上に表示されることになる。
【００７５】
たとえば、位相勾配分布画像データ生成処理部４１は、位相勾配データの値が正であれば第１の色データ（たとえば青色データ）を位相勾配分布画像データとし、位相勾配データの値が負であれば第１の色データとは異なる第２の色データ（たとえば赤色データ）を位相勾配分布画像データとするものであってもよい。この場合、入力装置１１から、子宮頚部から子宮体部に向かう方向（またはその正反対方向）を位相勾配算出方向として指示しておけば、頚部から体部に向かう蠕動と、体部から頚部に向かう蠕動とを、色分けして表示することができる。
【００７６】
また、位相勾配分布画像データ生成処理部４１は、位相勾配データの値が正であれば第１の色データ（たとえば青色データ）を割り当て、位相勾配データの値が負であれば第１の色データとは異なる第２の色データ（たとえば赤色データ）を割り当てるとともに、さらに、位相勾配データの絶対値を輝度値（画素値）とする位相勾配分布画像データを生成するものであってもよい。このようにすれば、各部の蠕動の方向を色分けして表すとともに、各部の蠕動の速さを輝度によって表すことができる。
【００７７】
画像オーバレイ処理部４２は、原子宮ＭＲ画像データ記憶部５１から子宮ＭＲ画像データ（たとえば、アフィン変換処理の際の基準とした１枚の画像）を読み出し、位相勾配分布画像データ記憶部６０から位相勾配分布画像データを読み出し、これらをオーバレイしたオーバレイ画像データを生成して、オーバレイ画像データ記憶部６１に格納する。すなわち、画像オーバレイ処理部４２は、子宮ＭＲ画像データにおいて、トリミング処理部３２によって切り出された矩形領域に対応する領域に、位相勾配分布画像データを重ね合わせる。画像表示制御部７０は、オーバレイ画像データ記憶部６１からオーバレイ画像データを読み出し、これをディスプレイ１３に画像表示する。これによって、子宮ＭＲ画像（原画像）に位相勾配分布画像が重ね合わせられて表示されるので、子宮の蠕動の様子を、その位置とともに一目で把握できるようになる。
【００７８】
むろん、画像表示制御部７０により、位相勾配分布画像データ記憶部６０に記憶された画像データをディスプレイ１３に画像表示させることもできる。このときには、位相勾配分布画像のみがディスプレイ１３に表示される。
速度測定点入力受付け部４３は、位相画像、位相勾配分布画像またはオーバレイ画像がディスプレイ１３に表示されている状態で、入力装置１１からの画像中の２点の入力を受け付ける。この入力は、当該２点間の速度を計測すべき速度測定点として受け付けられる。
【００７９】
位相差演算部４４は、速度測定点入力受付け部４３によって受け付けられた２点の画素の位相データを位相画像データ記憶部５８から読み出し、それらの差を位相差データとして演算して、位相差データ記憶部６２に格納する。
蠕動速度演算部４５は、速度測定点入力受付け部４３によって受け付けられた２つの速度測定点間の距離と、位相差データ記憶部６２に記憶された位相差データとに基づき、前記２点間の速度を蠕動速度として求め、この蠕動速度を表す蠕動速度データを蠕動速度データ記憶部６３に格納する。
【００８０】
画像表示制御部７０は、蠕動速度データ記憶部６３に格納された蠕動速度データをディスプレイ１３上に表示（たとえば数値表示）する。
図７は、子宮ＭＲ画像データを処理して子宮の蠕動を可視化するための処理の流れをまとめて示すフローチャートである。ＭＲＩによって得られた複数枚の子宮ＭＲ画像データに対して、画像間の子宮の位置および形状の変動を補償するための幾何学変換であるアフィン変換処理が行われ（ステップＳ１）、さらに、子宮部分の画像を切り出すためのトリミングが行われる（ステップＳ２）。次に、トリミングされた複数枚の子宮ＭＲ画像データを対象画像データとして以下の処理が行われる。
【００８１】
すなわち、複数枚の子宮ＭＲ画像データに対して、時間軸方向のフーリエ変換処理（ＦＦＴ）が実行され、個々の画素について、複数の周波数成分に関する振幅データおよび位相データが生成される（ステップＳ３）。次いで、個々の画素に関して、最大の振幅データを有する周波数成分（最大振幅周波数成分）が抽出され（ステップＳ４）、この最大振幅周波数成分に対応した色データを画素値とした疑似カラー画像である最大振幅周波数分布画像データが作成され、これに対応する画像（最大振幅周波数分布画像）がディスプレイ１３に表示される（ステップＳ５）。
【００８２】
操作者（典型的には、画像診断を行う専門医）は、ディスプレイ１３に表示された最大振幅周波数分布画像に基づき、子宮の蠕動に対応すると考えられる周波数成分を選択して、選択指示を入力装置１１から入力する（ステップＳ６）。
これにより、選択された周波数成分の振幅画像データおよび位相画像データが作成され、これらにそれぞれ対応する振幅画像および位相画像がディスプレイ１３に表示される（ステップＳ７，Ｓ８）。操作者は、表示された振幅画像および位相画像から、選択した周波数成分が子宮蠕動に対応しているかどうかを判断する（ステップＳ９）。当該選択した周波数が子宮蠕動を反映していると考えられる場合には、操作者は、以降の処理の実行を、入力装置１１から指示する（ステップＳ１０）。これにより、ステップＳ１１以降の処理が実行される。一方、前記選択された周波数成分が、子宮蠕動を反映していないと考えられる場合には、操作者は、別の周波数成分の選択指示を入力装置１１から入力する（ステップＳ６）。これにより、新たに選択された当該別の周波数成分に関して、振幅画像データおよび位相画像データが作成され、対応する振幅画像および位相画像がディスプレイ１３にそれぞれ表示される（ステップＳ７，Ｓ８）。
【００８３】
ステップＳ１１では、操作者は、子宮蠕動を反映していると判断された周波数成分の位相画像上の２点を入力装置１１から指定することにより、位相勾配の算出方向を指定する。これにより、当該方向に沿って、位相画像中の画素間の輝度勾配が位相勾配データとして求められる（ステップＳ１２）。さらに、この位相勾配データの符号に対応した２色の色データを画素値とする位相勾配分布画像データが作成され、対応する位相勾配分布画像がディスプレイ１３に表示される（ステップＳ１３）。さらに、必要に応じて、使用者は、子宮ＭＲ画像データと位相画像分布画像とのオーバレイ画像の表示を、入力装置１１から指示する（ステップＳ１４）。これに応答して、両画像をオーバレイしたオーバレイ画像データが作成されて、対応するオーバレイ画像がディスプレイ１３上に表示される（ステップＳ１５）。
【００８４】
子宮蠕動の速度を求める場合には、位相画像（ステップＳ８）、位相勾配分布画像（ステップＳ１３）およびオーバレイ画像（ステップＳ１５）のいずれかにおいて、操作者は、入力装置１１から画像中の２点を速度測定点として入力する（ステップＳ１６）。これに応答して、入力された２つの速度測定点間の位相データの差（位相差）が求められ（ステップＳ１７）、この位相差を用いて、蠕動運動の速度が求められ（ステップＳ１８）、そして、求められた蠕動運動の速度がディスプレイ１３上に表示される（ステップＳ１９）。
［シミュレーション実験］
子宮ＭＲ画像に人工的に内膜運動を付け加えた画像を作成し、その画像に対して前述の画像処理を行った。
【００８５】
図８は、シミュレーション実験のために作成した６４枚の時系列子宮ＭＲ画像を説明するための画像例である。１枚の対象画像（子宮ＭＲ画像）上の位置Ａに黒い点を付け加えて１枚目の子宮ＭＲ画像とし、２枚目以降の画像には同じ対象画像上に、順次、位置Ｂに向かって位置をずらせながら同じく黒い点を付け加え、８枚目で黒い点が位置Ｂに到達するようにした。９枚目以降の画像についても同様に、同じ対象画像を用い、８枚で位置Ａから位置Ｂまで移動するように黒い点を画像中に付け加えた。こうして、８枚で１周期となる模擬蠕動運動を付け加えた６４枚の子宮ＭＲ画像データを作成した。
【００８６】
位置ＡおよびＢの画素の輝度値の時間変化を図９に示す。横軸は、時間であるが、１枚目から６４枚目までの画像の番号で表されている。この図より、黒い点が移動することによって、位置Ａ，Ｂの画素の輝度値が変化している様子が観測できる。黒い点が位置Ａ，Ｂに位置する時間が異なるので、位置Ａ，Ｂの画素にそれぞれ対応した２つの曲線（折れ線）が、時間軸方向にずれている。
【００８７】
図８の例において、位置Ａの座標は（１１９，８７）、位置Ｂの座標は（９１，１１１）であった。この場合、２点間の距離は４４．６５（ｍｍ）となる。つまり、（８−１）×２（ｓｅｃ）で４４．６５（ｍｍ）を移動するため、ここでシュミレーションデータとして作成した運動の速度は、４４．６５（ｍｍ）／１４（ｓｅｃ）＝３．１９（ｍｍ／ｓｅｃ）となる。
【００８８】
前述のような６４枚の子宮ＭＲ画像データに対して時間軸方向にフーリエ変換を施し、各画素の最大振幅の周波数成分を求め、最大振幅周波数分布画像を作成した。これを図１０（(a)はカラー画像、(b)は(a)のカラー画像の第８周波数の色部分を強調した画像）に示す。この最大周波数分布画像から、第８周波数成分（＝６４／８）が運動を表現していることが分かる。この第８周波数成分に対する振幅画像および位相画像を図１１および図１２にそれぞれ示す。
【００８９】
振幅画像（図１１）において、高輝度となっている場所が、人工画像（黒い点を付加した６４枚の子宮ＭＲ画像）で黒い点が移動した場所である。位相画像（図１２）では、所定のしきい値未満の振幅の画素の画素値は「０」に固定されるので、図１１の振幅画像において高輝度となっている領域と同じ領域に有意な画素値を持つ画素が現れている。この位相画像の画素値（位相の値）を用いることで、２点間の運動の速度（向きと速さ）を解析できる。つまり、この場合、高輝度な位置Ａの方から低輝度な位置Ｂに向かう運動として検出することができる。
【００９０】
位相画像中の位置Ａと位置Ｂは、人工画像を作成した際の運動の始点と終点である。この２点の距離と位相差より、運動の速度を求める。位置（１１９，８７）と位置Ｂ（９１，１０８）の間の距離は４４．６５（ｍｍ）であり、位相はそれぞれ３．１４（ｒａｄ），−２．３６（ｒａｄ）であった。したがって、位置Ａ，Ｂの画素間の位相差は５．５０（ｒａｄ）である。よって、運動の速度は、３．１９（ｍｍ／ｓｅｃ）となり、設定した速度に等しい結果が得られた。
【００９１】
位相勾配を求める方向を、位置Ａから位置Ｂに向かう方向に指定して、各画素ごとに隣接画素との輝度値の差（位相勾配）を求め、正の位相勾配（位置Ａから位置Ｂに向かって位相が増加する場合に対応）に対して青色の色データを対応付け、負の位相勾配（位置Ｂから位置Ａに向かって位相が増加する場合に対応）に対して赤色の色データを対応付け、これらの色データを画素値とする位相勾配分布画像を作成すると、図１３に示す画像が得られる。すなわち、位置Ａ，Ｂ間の領域が青色に着色された画像が得られる。
［実験１（排卵期の画像）］
図１４に示す対象画像（６４枚の子宮ＭＲ画像。図１４にはその中の代表的な１枚を示す。）は或る被験者の排卵期の子宮をＭＲＩによって撮影したものである。この対象画像に本実施形態の画像処理を適用した。
【００９２】
図１５(a)に最大振幅周波数分布画像（カラー画像）を示す。図１５(b)には図１５(a)のカラー画像中の第３周波数の色部分を強調した画像を示し、図１５(c)には図１５(a)のカラー画像中の第６周波数の色部分を強調した画像を示し、図１５(d)には図１５(a)のカラー画像中の第７周波数の色部分を強調した画像を示す。
図１５(a)の画像から、子宮内膜の運動の周波数が一定でないこと、すなわち、内膜の運動にはいくつかの異なる周波数の運動があることが分かる。この例では、最大振幅周波数分布画像から一意に適切な周波数成分（内膜の蠕動に対応する周波数成分）を特定できない。そこで、振幅画像と位相画像をいくつかの周波数成分について作成した。その結果、第７周波数成分の振幅画像および位相画像がそれぞれ図１６および図１７に示すとおりとなり、子宮内膜の形状の画像部分が形成された。したがって、第７周波数成分が子宮蠕動に対応するものであると判断した。
【００９３】
図１７において、内膜部分が、子宮の頚部から体部に向かって徐々に暗くなっているため、対象画像子宮の内膜の運動方向は、頚部から体部に向かう方向であることが分かる。
図１７の４つの点ａ，ｂ，ｃ，ｄの輝度値の時間軸方向の変化を図１８および図１９に示す。図１８は原データ（フーリエ変換前のデータ）、図１９は直流成分と第７周波数成分に逆フーリエ変換を適用したものである。図１９を見ると、４つの点ａ，ｂ，ｃ，ｄの輝度値変化が時間軸方向にずれていることが分かる。この時間軸方向のずれが、図１７の位相画像に輝度値の変化として現れている。
【００９４】
排卵期には、精子の輸送のため、子宮の内膜が頚部から体部へ運動していると考えられており、排卵期の子宮の画像を扱った本実験では、そのことが確かめられた。
内膜上部の速度は図１７中の点ａと点ｄの位相から求めた。点ａの座標は（１１４，７８）で、点ｄの座標は（７５，９１）であったので、２点間の距離は５０．２１（ｍｍ）となる。位相はそれぞれ１．８９（ｒａｄ）、−１．５７（ｒａｄ）であったので、位相差は３．４６（ｒａｄ）である、よって、速度は５．３２（ｍｍ／ｓｅｃ）となる。内膜下部の速度は、図１７の画像中の点６と点７から求めた。点６の座標は（１０９，９４）で、点７の座標は（８９，１０９）であったので、２点間の距離は３１．８９（ｍｍ）となる。位相はそれぞれ１．０９（ｒａｄ），−１．７１（ｒａｄ）だったので、位相差は２．８０（ｒａｄ）である。よって、速度は４．１７（ｍｍ／ｓｅｃ）となる。
【００９５】
本手法により解析された運動が、画像列を連続で表示したときに視覚で捉えられる運動と同じ運動を捉えているかどうかを評価するために、画像中で観察できる運動の速度および位相から運動の速度を求めた。
図２０は、前記の実験の対象画像の画像列の２８枚目の子宮ＭＲ画像を示し、図２１は当該画像列の３７枚目の子宮ＭＲ画像を示す。この実験の対象画像では、図２０に示す２８枚目から図２１に示す３７枚目の画像の間で、黒い影（図中の矢印で示す部分）が、内
膜上部では点Ｃ（１２３，７３）から点Ｄ（５７，１０３）まで８９．６５（ｍｍ）移動し、内膜下部では点Ｅ（１１９，８７）から点Ｆ（７１，１１５）まで６９．６７ｍｍ移動していた。経過時間は１８（ｓｅｃ）（＝２sec×９）であるため、内膜上部と下部の運動の速度は、それぞれ４．９８（ｍｍ／ｓｅｃ）、３．８７（ｍｍ／ｓｅｃ）であった。 下記表１に示すように、本手法により求めた速度と、画像中で捉えられる運動の速度がそれぞれほぼ近い値となっており、本手法で解析した結果と、目視で検出した結果とが、同等な運動を捉えていると考えられる。
【００９６】
【表１】
【００９７】
よって、本手法が、一般には手作業で計測することが困難な子宮の蠕動運動を、画像解析によって求めることができる有効な手法であることが確認された。
位相勾配を求める方向を、図１７の点ａから点ｄ向かう方向に指定して、各画素ごとに隣接画素との輝度値の差（位相勾配）を求め、正の位相勾配（点ａから点ｄに向かって位相が増加する場合に対応）に対して青色の色データを対応付け、負の位相勾配（点ｄから点ａに向かって位相が増加する場合に対応）に対して赤色の色データを対応付け、これらの色データを画素値とする位相勾配分布画像を作成すると、図２２（(a)はカラー画像、(b)は(a)のカラー画像の青色部分を強調した図。右下部の矢印は、左下に向かうものが青色矢印であり、右上に向かうものが赤色矢印である。）に示す画像が得られる。すなわち、子宮内膜に対応する部分が青色に着色された画像が得られ、子宮頚部から体部に向かう蠕動運動が明瞭に可視化されている。
【００９８】
このカラー画像を原データに重ね合わせることにより、図２３に示すオーバレイ画像が得られる（(a)はカラー画像、(b)は(a)のカラー画像の青色領域を白抜きで強調表示した画像）。
［実験２（月経期の画像］
図２４に示す対象画像（６４枚の子宮ＭＲ画像。図２４にはそのなかの代表的な１枚を示す。）は、実験１の場合と同じ被験者の月経期の画像である。この画像に前記実施形態の手法を適用し、排卵期と月経期の運動の違いを検証した。
【００９９】
図２５(a)に最大振幅周波数分布画像を示す。参考のため、図２５(b)には図２５(a)のカラー画像中の第２周波数の色部分を強調した画像を示し、図２５(c)には図２５(a)のカラー画像中の第３周波数の色部分を強調した画像を示す。
図２５(a)の画像から、主な運動の周波数は第２周波数成分であると判断した。この第２周波数成分の振幅データおよび位相データを用いて作成した振幅画像および位相画像を図２６および図２７にそれぞれ示す。
【０１００】
図２７において、内膜部分が、子宮の体部から頚部に向かって徐々に暗くなっていることが観察できる。したがって、運動方向は、排卵期とは逆で、体部から頚部に向かう方向であると判断できる。月経期の子宮は、月経血の排出のため、体部から頚部に向かって運動していると考えられており、この実験では、そのことが確かめられた。
内膜上部の運動の速度を図２７中の点９と点８の画素の位相から求め、内膜下部の運動の速度を図２７中の点ｈと点ｅの位相から求めた。求めた速度を表２に示す。
【０１０１】
【表２】
【０１０２】
位相勾配を求める方向を、図２７の点８から点９向かう方向に指定して、各画素ごとに隣接画素との輝度値の差（位相勾配）を求め、正の位相勾配（点８から点９に向かって位相が増加する場合に対応）に対して青色の色データを対応付け、負の位相勾配（点９から点８に向かって位相が増加する場合に対応）に対して赤色の色データを対応付け、これらの色データを画素値とする位相勾配分布画像を作成すると、図２８(a)に示すカラー画像が得られる（(b)は(a)のカラー画像中の赤色部分を強調表示した画像。右下部の矢印は、左下に向かうものが青色矢印であり、右上に向かうものが赤色矢印である。）。すなわち、子宮内膜に対応する部分が赤色に着色された画像が得られ、子宮体部から頚部に向かう蠕動運動が明瞭に可視化されている。
【０１０３】
このカラー画像を原データに重ね合わせることにより、図２９に示すオーバレイ画像が得られる（ (a)はカラー画像、(b)は(a)のカラー画像の赤色領域を白抜きで強調表示した画像）。
［実験３（黄体期の画像］
図３０に示す対象画像（６４枚の子宮ＭＲ画像。図３０には代表的な１枚のみを示す。）は、或る被験者の黄体期の子宮の画像である。この画像の内膜の運動を解析し、黄体期の子宮内膜が、排卵期、月経期とはどのように異なるかを考察した。
【０１０４】
図３１(a)に、当該対象画像の最大振幅周波数分布画像を示す。参考のため、図３１(b)には図３１(a)のカラー画像中の第１周波数の色部分を強調した画像を示し、図３１(c)には図３１(a)のカラー画像中の第２周波数の色部分を強調した画像を示す。
図３１(a)の画像から、主要な運動の周波数成分（各部で最大振幅をとる周波数成分）は、第１周波数成分であると判断し、その振幅データおよび位相データに基づいて、振幅画像（図３２）および位相画像（図３３）をそれぞれ作成した。
【０１０５】
図３２の振幅画像には、内膜部分に高輝度な場所が明確には存在していない。また、図３３から、内膜部分の輝度値が一定であるため、運動による位相の変化は存在していないことが分かる。
図３４は、図３０中の４点ｑ，ｒ，ｓ，ｔの原データ（輝度値の時間変化）であり、図３５は、フーリエ変換によって得られた点ｑ，ｒ，ｓ，ｔの周波数スペクトル（各周波数成分の振幅データ）である。図３４から、時間軸方向の周期的な輝度値の変化がないことがわかる。また、図３５の振幅データから、最大振幅をとる周波数成分が明確ではなく、どの周波数成分の振幅の値も小さいことが分かる。
【０１０６】
これらの結果より、黄体期では、排卵期や月経期のような内膜運動は行われていないと判断できる。黄体期(分泌期) には、妊卵の着床や初期の妊卵の保持のため、内膜は蠕動運動を行っていないと考えられており、本実験による解析結果と一致する。
位相勾配を求める方向を、対称画像の右上角から左下角に向かう方向に指定して、各画素ごとに隣接画素との輝度値の差（位相勾配）を求め、正の位相勾配（右上角から左下角に向かって位相が増加する場合に対応）に対して青色の色データを対応付け、負の位相勾配（左下角から右上角に向かって位相が増加する場合に対応）に対して赤色の色データを対応付け、これらの色データを画素値とする位相勾配分布画像を作成すると、図３６に示す画像が得られる。この画像には、蠕動運動が現れていない。
［実験４（筋層の運動の解析）］
この実験では、子宮の内膜の外側を覆うように存在する子宮の筋層の運動に注目した。筋層は厚い筋肉から成り立っており、内膜と同様に周期的な運動をしていると考えられている。
【０１０７】
対象画像は図１４の場合（実験１）と同じである。したがって、最大振幅周波数分布画像は図１５(a)のとおりである。
この図１５(a)の最大振幅周波数分布画像からは、運動の周波が一定ではないと認められるため、筋層の運動を表す周波数成分を一意に決定することができない。そこで、いくつか周波数成分について、振幅画像と位相画像を作成し、図３７および図３８に振幅画像および位相画像をそれぞれ示す第６周波数成分が筋層に対応しているものと判断した（前述のとおり、第７周波数成分は内膜の蠕動に対応している）。
【０１０８】
図３８の位相画像より、運動方向は頚部から体部に向かう方向であると判断できる。筋層の運動は、内膜の運動よりも、画像中から手作業で捉えるのは一層困難であるが、本手法により解析が可能であることが確認できた。
排卵期の内膜の運動は頚部から体部に向かう方向の運動であったが、この実験により、筋層の運動方向もから体部に向かう方向であるという結果が得られた。つまり、内膜と筋層の運動方向は同じであるという結果が得られた。
【０１０９】
運動の速度を、図３８中の点ｕと点ｘから求めた。結果を表３に示す。
【０１１０】
【表３】
【０１１１】
次に、画像中で捉えられる筋層の運動から、運動の速度を手作業で求めた。画像列の１８枚目（図３９）から２４枚目（図４０）の間で、黒い影が、筋層を中の点Ｓから点Ｔまで移動するのが認められた。この運動から目視で求めた速度を表４に示す。
【０１１２】
【表４】
【０１１３】
画像中から目視で捉えられる運動の速度と、本手法により捉えた運動の速度がほぼ同じ値となっているため、筋層の蠕動の解析にも本手法が有効であると言える．
［他の実施形態］
前述の実施形態では、子宮の蠕動に対応した周波数成分の選択のために、最大振幅周波数分布画像がディスプレイ１３に表示され、操作者の判断によって、振幅画像および位相画像を作成すべき周波数成分が選択されているが、このような周波数成分の選択を自動化してもよい。具体的には、たとえば、最大振幅周波数分布画像中における個々の周波数成分の画素の数（すなわち、面積）を計算し、この画素数（面積）が上位の（たとえば最大の）周波数成分を自動で選択する周波数成分選択処理部４６（図６参照）を設けるようにしてもよい（面積順位基準選択手段）。また、周波数成分選択処理部４６は、所定のしきい値以上の最大振幅データに対応した１つまたは複数個の周波数成分を、子宮の蠕動に対応する蓋然性が高い周波数成分として自動選択する処理を行うものであってもよい。
【０１１４】
また、前述の実施形態では、算出すべき位相勾配の方向を入力装置１１から操作者が指定する構成としているが、位相勾配の方向は、予め定めた一定の方向（右上角から左下角に向かう方向、右から左に向かう方向など）であってもよい。
その他、特許請求の範囲に記載された事項の範囲で種々の設計変更を施すことが可能である。
【図面の簡単な説明】
【０１１５】
【図１】子宮ＭＲ画像（矢状断面）の一例を示す。
【図２】子宮ＭＲ画像の例を示す。
【図３】時間軸方向の輝度値の変化をプロットした例を示す。
【図４（ａ）】最大振幅周波数分布画像（カラー画像）の一例を示す。
【図４（ｂ）】図４(a)のカラー画像の第１周波数の色部分を強調表示した画像を示す。
【図４（ｃ）】図４(a)のカラー画像の第２周波数の色部分を強調表示した画像を示す。
【図４（ｄ）】図４(a)のカラー画像の第６周波数の色部分を強調表示した画像を示す。
【図５】この発明の一実施形態に係る画像診断支援システムの基本的なハードウェア構成を説明するためのブロック図である。
【図６】前記画像処理装置として機能するコンピュータ本体の機能的な構成を説明するためのブロック図である。
【図７】子宮ＭＲ画像データを処理して子宮の蠕動を可視化するための処理の流れをまとめて示すフローチャートである。
【図８】シミュレーション実験のために作成した６４枚の時系列子宮ＭＲ画像を説明するための画像例である。
【図９】図８の位置ＡおよびＢの画素の輝度値の時間変化を示す。
【図１０】図８の子宮ＭＲ画像についての最大振幅周波数分布画像を示す。
【図１１】振幅画像を示す。
【図１２】位相画像を示す。
【図１３】位相勾配分布画像を示す。
【図１４】排卵期の子宮ＭＲ画像の例を示す。
【図１５（ａ）】図１４の画像の最大振幅周波数分布画像を示す。
【図１５（ｂ）】図１５(a)のカラー画像の第３周波数の色部分を強調表示した画像を示す。
【図１５（ｃ）】図１５(a)のカラー画像の第６周波数の色部分を強調表示した画像を示す。
【図１５（ｄ）】図１５(a)のカラー画像の第７周波数の色部分を強調表示した画像を示す。
【図１６】子宮内膜の蠕動に対応した振幅画像を示す。
【図１７】子宮内膜の蠕動に対応した位相画像を示す
【図１８】図１７の４つの点ａ，ｂ，ｃ，ｄの輝度値の時間軸方向の変化を示す（原データ）。
【図１９】図１７の４つの点ａ，ｂ，ｃ，ｄの輝度値の時間軸方向の変化を示す（直流成分と第７周波数成分に逆フーリエ変換を適用したもの）。
【図２０】図１４の子宮ＭＲ画像列の２８枚目の子宮ＭＲ画像を示す。
【図２１】図１４の子宮ＭＲ画像列の３７枚目の子宮ＭＲ画像を示す。
【図２２】位相勾配分布画像を示す。
【図２３】オーバレイ画像を示す。
【図２４】月経期の子宮ＭＲ画像の例を示す。
【図２５】最大振幅周波数分布画像を示す。
【図２６】子宮内膜の蠕動に対応した振幅画像を示す。
【図２７】子宮内膜の蠕動に対応した位相画像を示す。
【図２８】位相勾配分布画像を示す。
【図２９】オーバレイ画像を示す。
【図３０】黄体期の子宮ＭＲ画像の例を示す。
【図３１】最大振幅周波数分布画像を示す。
【図３２】振幅画像を示す。
【図３３】位相画像を示す。
【図３４】図３０中の４点ｑ，ｒ，ｓ，ｔの原データ（輝度値の時間変化）を示す。
【図３５】図３０中の４点ｑ，ｒ，ｓ，ｔの周波数スペクトル（各周波数成分の振幅データ）である。
【図３６】位相勾配分布画像を示す。
【図３７】子宮筋層の蠕動に対応する振幅画像を示す。
【図３８】子宮筋層の蠕動に対応する位相画像を示す。
【図３９】図１４の子宮ＭＲ画像列の１８枚目の画像を示す。
【図４０】図１４の子宮ＭＲ画像列の２４枚目を示す。
【符号の説明】
【０１１６】
１０ コンピュータ本体
１１ 入力装置
１１Ａ キーボード
１１Ｂ マウス
１３ ディスプレイ
１４ プリンタ
１５ ＣＰＵ
１６ ＲＯＭ
１７ ＲＡＭ
１８ ハードディスクドライブ
１９ ＣＤ−ＲＯＭドライブ
２０ 入力制御部
２１ ディスプレイインタフェース
２２ プリンタインタフェース
２３ バス
３０ 画像処理演算部
３１ アフィン変換処理部
３２ トリミング処理部
３３ フーリエ変換処理部
３４ 最大振幅周波数成分抽出処理部
３５ 最大振幅周波数分布画像データ生成処理部
３６ 周波数成分選択入力受付け部
３７ 振幅画像データ生成処理部
３８ 位相画像データ生成処理部
３９ 位相勾配算出方向入力受付け部
４０ 位相勾配データ生成処理部
４１ 位相勾配分布画像データ生成処理部
４２ 画像オーバレイ処理部
４３ 速度測定点入力受付け部
４４ 位相差演算部
４５ 蠕動速度演算部
４６ 周波数成分選択処理部
５０ データ記憶部
５１ 原子宮ＭＲ画像データ記憶部
５２ アフィン変換後子宮ＭＲ画像データ記憶部
５３ 対象子宮ＭＲ画像データ記憶部
５４ フーリエ変換データ記憶部
５５ 最大振幅周波数成分データ記憶部
５６ 最大振幅周波数分布画像データ記憶部
５７ 振幅画像データ記憶部
５８ 位相画像データ記憶部
５９ 位相勾配データ記憶部
６０ 位相勾配分布画像データ記憶部
６１ オーバレイ画像データ記憶部
６２ 位相差データ記憶部
６３ 蠕動速度データ記憶部
７０ 画像表示制御部
【技術分野】
【０００１】
この発明は、時系列に従って子宮を撮影して得た複数枚の画像データに基づいて子宮の蠕動運動を解析して可視化するための画像処理装置、画像処理方法および画像診断支援システム、ならびにそのためのコンピュータプログラムおよびそれを記録した記録媒体に関する。
【背景技術】
【０００２】
婦人科領域の画像診断では、超音波画像、Ｘ線撮影画像およびＭＲ（Magnetic Resonance; 磁気共鳴) 画像が用いられている。超音波画像には死角が存在し、分解能が低く、組織特異性が描画されないという欠点があり、また、Ｘ線撮影にはＸ線被曝の問題がある。これに対して、ＭＲ画像を用いた診断は、Ｘ線を使用しないため、妊娠中の患者に対しても安全に撮影が行え、かつ、優れた組織コントラストにより、機能的変化および形態的変化を顕著に描画することができる。
【０００３】
従来型のＭＲＩ（Magnetic Resonance Imaging）は、撮影時間が長く、子宮内膜が波打つような反復した収縮運動や子宮内膜の外側の筋層の反復運動のような子宮蠕動運動を捉えることは不可能であった。そのため、画像診断において、子宮の蠕動は全く考慮されず、子宮はあたかも静止した臓器であるかのように扱われてきた。
しかし、近年になって、ＨＡＳＴＥ（half-Fourier-acquisition single-shot turbo spin-echo）等の高速撮影法が開発されるに至り、ＭＲ画像の撮影に要する時間が短縮され、これに伴い、子宮蠕動運動を直接描画する動態撮像が可能となってきた（非特許文献１、２）。
【０００４】
女性生殖器は、年齢、月経周期およびホルモン等の影響を受ける。そのため、子宮は、他の臓器とは異なり、正常時でもその形状は一定ではない。すなわち、子宮は、排卵期、月経期、黄体期に、それぞれ異なる性質を持った運動を行っていると考えられている。そこで、子宮の蠕動運動を解析することにより、月経、妊能に関する各患者の症状に適した治療法の発見につながることが期待されている。
【非特許文献１】植田隆史他著「骨盤部（女性）撮像に必要な解剖・生理機能」、INNERVISION(18・9)、２００３年
【非特許文献２】富樫かおり著「ＭＲＩによる婦人科画像診断」、日本放射線技術学会雑誌第５９巻第８号、２００３年８月、ｐｐ８９５−９０３
【発明の開示】
【発明が解決しようとする課題】
【０００５】
しかし、現在のＭＲ画像による子宮蠕動運動の分析は、生のＭＲ画像の肉眼による目視観察に依存している。ところが、生のＭＲ画像中に現れる子宮の蠕動は、子宮内膜に現れる微少な波状突起の移動や、子宮筋層部分内における微妙な陰影（低輝度部）の移動として観察されるに過ぎない。より具体的に説明すると、骨盤部を撮影したＭＲ画像中では、周囲の臓器の活発な運動が現れており、このような周囲の臓器の運動の影響を受けて子宮の形状も絶えず変動している。このような状況の画像中では、子宮の蠕動は、ノイズに近い微弱な運動として観察されるに過ぎず、熟練した専門医師による極めて注意深い観察を必要とする。
【０００６】
そのため、ＭＲ画像に基づく子宮蠕動運動の評価は主観的にならざるをえず、適正な診断には、専門医師であっても、相当の熟練を要するという問題がある。
そこで、この発明の目的は、子宮の蠕動運動を分かりやすく可視化することができ、これにより、子宮の画像診断を支援することができる画像処理装置、画像処理方法および画像診断支援システム、ならびにそのためのコンピュータプログラムおよびそれを記録した媒体を提供することである。
【課題を解決するための手段】
【０００７】
上記の目的を達成するための請求項１記載の発明は、子宮の蠕動運動を解析して可視化するための画像処理装置であって、子宮の蠕動の周期よりも短い周期（たとえば、２〜３秒間隔）で時系列に従って撮影した複数枚の子宮ＭＲ（磁気共鳴）画像に対応した複数枚の子宮ＭＲ画像データ（個々の画素の輝度データで構成されるもの。とくに、いわゆる矢状断面画像）を記憶するための画像データ記憶手段と、この画像データ記憶手段に記憶された複数枚の子宮ＭＲ画像データに対して時間軸方向のフーリエ変換を行い、子宮ＭＲ画像中の各画素について複数の周波数成分（輝度の時間変化の周波数成分）の位相データを生成するフーリエ変換手段と、このフーリエ変換手段によって生成された少なくとも１つの周波数成分の位相データを用いて、子宮ＭＲ画像各部における当該周波数成分の所定方向に関する位相勾配（当該所定方向に隣接する画素との位相データの差分。微分値）を算出し、その位相勾配を表す位相勾配データを生成する位相勾配データ生成手段と、この位相勾配データ生成手段によって生成された位相勾配データに基づき、画像各部を位相勾配に応じて着色した位相勾配分布画像を表す位相勾配分布画像データを、子宮の蠕動を表す画像データ（蠕動の方向および／または速さを表すもの）として生成する位相勾配分布画像データ生成手段とを含むことを特徴とする子宮蠕動可視化のための画像処理装置である。
【０００８】
子宮の蠕動は、子宮内膜における波状の微少突起のほぼ一定速度での周期的な移動と、子宮筋層における局部的な状態変化部位のほぼ一定速度での周期的な移動とに大別される。一般に、排卵期には、子宮内膜および子宮筋層における蠕動方向は、精子の輸送のために、子宮の頚部から体部へと向かう方向である。月経期には、子宮内膜および子宮筋層における蠕動方向は、月経血の排出のために、子宮体部から子宮頚部へと向かう方向である。さらに、黄体期（分泌期）では、妊卵の着床や初期の妊卵の保持のため、内膜は蠕動を行っていないと考えられている。
【０００９】
この発明では、時系列に従う複数枚の子宮ＭＲ画像データに対して時間軸方向のフーリエ変換が施される。これにより、子宮ＭＲ画像を構成する個々の画素について、画素値（通常は輝度値）の時間変化が複数の周波数成分に分解され、個々の周波数成分の位相データが生成される。前述のように、子宮の蠕動は、ほぼ一定速度での子宮内膜および子宮筋層の周期的な運動であるので、この蠕動に対応した周波数成分に着目することによって、子宮の蠕動を解析することができる。
【００１０】
そこで、この発明では、フーリエ変換によって生成された少なくとも１つの周波数成分の位相データが用いられ、子宮ＭＲ画像各部における当該周波数成分の所定方向に関する位相勾配を表す位相勾配データが算出される。子宮の蠕動に対応した周波数成分の位相データに関して所定方向の勾配（位相勾配）を算出すると、この位相勾配は、子宮各部の運動速度（蠕動の速さおよび方向）を表す。このような位相勾配データに応じた色データを画素値とする位相勾配分布画像を形成すると、子宮各部の運動（方向および／または速さ）を一目で把握することができる。この発明では、このような位相勾配分布画像を表す位相勾配分布画像データが生成される。
【００１１】
したがって、子宮の蠕動に対応した周波数成分の位相データに対応した位相勾配分布画像データを用いて表示装置に位相勾配分布画像を表示することにより、子宮蠕動を明瞭に可視化した画像を提供することができる。これにより、子宮に関する画像診断を効果的に支援することができる。
位相勾配を求める方向は、操作者が指定することとしてもよいし、所定の一定方向に予め設定しておいてもよい。
【００１２】
前記画像処理装置は、前記フーリエ変換手段によって生成された少なくとも１つの周波数成分の位相データを用いて、子宮ＭＲ画像各部の当該周波数成分の位相分布を表す位相分布データ（位相画像データ）を生成する位相分布データ生成手段（位相画像データ生成手段）を含んでいてもよい。この場合、前記位相勾配データ生成手段は、前記位相分布データ生成手段によって生成された位相分布データ（各画素の位相データ）に基づいて、子宮ＭＲ画像の各部の位相勾配を算出するものであってもよい。
【００１３】
請求項２記載の発明は、前記画像データ記憶手段に記憶された複数枚の子宮ＭＲ画像データに対して、子宮画像部分の位置および形状の変動を補償するための幾何学変換を施し、この幾何学変換処理後の子宮ＭＲ画像データを生成する幾何学変換手段をさらに含み、前記フーリエ変換手段は、前記幾何学変換手段による幾何学変換処理後の子宮ＭＲ画像データに対してフーリエ変換を実行するものであることを特徴とする請求項１記載の画像処理装置である。
【００１４】
ＭＲＩによって撮影される画像中で、子宮の位置および形状は時間経過とともに変動する。そこで、この発明では、複数枚の子宮ＭＲ画像データに対して、子宮画像部分の位置および形状を合わせるための幾何学変換を施し、その後に、フーリエ変換を行うようにしている。これにより、より確実に、子宮の蠕動に対応した画像の変化を抽出することができる。
【００１５】
請求項３記載の発明は、前記フーリエ変換手段は、前記位相データのほかに、子宮ＭＲ画像の各画素についての複数の周波数成分の振幅データを生成するものであり、このフーリエ変換手段が生成する振幅データに基づいて、少なくとも１つの周波数成分を選択する周波数成分選択手段をさらに含み、前記位相勾配データ生成手段は、前記周波数成分選択手段によって選択された周波数成分に関して子宮ＭＲ画像各部の位相勾配データを生成するものであることを特徴とする請求項１または２記載の画像処理装置である。
【００１６】
この構成では、フーリエ変換によって、子宮ＭＲ画像の各画素について、複数の周波数成分に関する位相データおよび振幅データが生成される。子宮ＭＲ画像中に見られる運動は、子宮の蠕動ばかりではなく、周辺の臓器からの影響による移動および変形などのノイズ的な運動も含まれている。フーリエ変換によって分解される複数の周波数成分のなかには、このようなノイズ的な運動を反映したものもある。
【００１７】
子宮を含む画像領域に注目すると、子宮蠕動に対応した周波数成分は、ノイズ的な運動に対応した周波数成分に比較すると、大きな振幅を持つ。そこで、フーリエ変換によって得られる個々の周波数成分の振幅データに注目することによって、子宮蠕動を表す周波数成分を選択することができる。そして、このようにして選択された周波数成分に関して子宮ＭＲ画像の各部の位相勾配データを求め、その分布を表す位相勾配分布画像データを作成すれば、この位相勾配分布画像データに対応する疑似カラー画像を表示装置に表示して、子宮の蠕動を分かりやすく可視化することができる。
【００１８】
なお、前記位相分布データ生成手段（位相画像データ生成手段）が備えられる場合には、この位相分布データ生成手段は、前記周波数成分選択手段によって選択された周波数成分に関する位相分布データを生成するものであることが好ましい。
請求項４記載の発明は、前記周波数成分選択手段は、子宮ＭＲ画像の各画素について、最大の振幅データに対応した最大振幅周波数成分を抽出する最大振幅周波数抽出手段と、この最大振幅周波数抽出手段によって個々の画素について抽出された最大振幅周波数成分の子宮ＭＲ画像内での分布を表す最大振幅周波数分布画像データ（たとえば、最大周波数成分に対応した色データを画素値とする画像データ）を生成する最大振幅周波数分布画像データ生成手段とを含むことを特徴とする請求項３記載の画像処理装置である。
【００１９】
この構成によれば、最大振幅周波数分布画像データに基づいて、表示装置に最大振幅周波数分布画像（たとえば、最大周波数成分の分布を表す疑似カラー画像）を表示すれば、子宮ＭＲ画像中での最大振幅周波数成分の分布を視覚的に把握することができる。この最大振幅周波数分布画像に基づいて、子宮の蠕動に対応した周波数成分を特定することができる。
【００２０】
より具体的には、請求項５に記載されているように、前記周波数成分選択手段は、さらに、操作者による周波数成分選択入力を受け付ける周波数成分選択入力受け付け手段を含むものであってもよい。
この画像処理装置を適用した画像診断支援システムでは、たとえば、操作者による入力操作が可能な入力装置と、前記最大振幅周波数分布画像データ生成手段が生成する最大振幅周波数分布画像データを画像表示する表示装置とが備えられる。操作者が、表示装置に表示された最大振幅周波数分布画像から判断して、入力装置からいずれかの特定の周波数成分の選択入力を行うことにより、この入力が前記周波数成分選択入力受け付け手段によって受け付けられ、その周波数成分が選択されることになる。
【００２１】
さらに具体的に説明すると、最大振幅周波数分布画像中において、子宮ＭＲ画像中の子宮の一部または全部に近似した形状を形成する周波数成分が存在すれば、このような周波数成分は、子宮蠕動に対応した周波数成分である可能性が高い。そこで、操作者は、このような周波数成分を選択して入力装置から入力することになる。
医療の現場における画像診断では、完全に自動化されたシステムよりも、操作者（一般的には専門医）による手動操作の余地を残すことが好まれる場合が多い。これは、診断の基礎となる画像は、個々の患者に依存する多様な画像であり、また、患者の状態によっても異なるため、画像の自動解析が必ずしも成功するとは限らず、必要に応じて専門医の判断に委ねるステップを設ける方が、信頼性の高い画像診断が可能であり、かつ、診断結果に対する安心感を得やすいからである。さらにまた、画像解析の一部のステップを操作者の判断に委ねることで、画像処理装置の演算負荷が著しく軽減されるという効果も得られる。
【００２２】
ただし、子宮蠕動に対応した周波数成分の選択の自動化が不可能であるというわけではなく、必要に応じて、周波数成分選択処理の自動化も可能である。
具体的には、請求項６に記載されているように、前記周波数成分選択手段は、さらに、前記最大振幅周波数分布画像データ生成手段が生成する最大振幅周波数分布画像データに基づき、当該最大振幅周波数分布画像データに対応した画像中に占める個々の周波数成分の面積の順位に基づいて、少なくとも１つの周波数成分を選択する面積順位基準選択手段を含むものであってもよい。
【００２３】
最大振幅周波数分布画像中では、子宮蠕動に対応した周波数成分の領域が支配的になる可能性が高い。したがって、最大振幅周波数分布画像中に占める面積が上位（たとえば第１位のみ。第１位〜第３位など）の周波数成分は、子宮蠕動に対応した周波数成分である可能性が高い。したがって、最大振幅周波数分布画像中での占有面積が上位の周波数成分を自動選択する前記の構成によって、操作者による判断を要することなく、子宮蠕動に対応した周波数成分を適切に選択することができる。
【００２４】
この構成は、フーリエ変換以降の処理対象となる子宮ＭＲ画像が、子宮近傍の領域を切り出した画像である場合にとくに効果的である。
また、請求項７に記載されているように、前記周波数成分選択手段が、所定のしきい値以上の最大振幅データに対応した周波数成分を選択するものである場合も、子宮蠕動に対応した周波数成分の自動選択が可能である。すなわち、子宮蠕動に対応した周波数成分はある程度の振幅を示すと考えられるので、所定のしきい値未満の最大振幅データに対応した周波数成分はノイズとみなせる。また、一定以上の振幅を持たない周波数成分については、その位相データに基づく解析自体があまり意味を持たなくなるから、無意味な画像解析処理を省くためにも、最大振幅データが前記しきい値に達しない周波数成分は、処理の対象外とすることが好ましい。
【００２５】
請求項８記載の発明は、前記フーリエ変換手段は、前記位相データのほかに、子宮ＭＲ画像の各画素についての複数の周波数成分の振幅データを生成するものであり、前記フーリエ変換手段によって生成された少なくとも１つの周波数成分の振幅データを用いて、子宮ＭＲ画像各部の当該周波数成分の振幅分布を表す振幅画像を表す振幅画像データを生成する振幅画像データ生成手段をさらに含むことを特徴とする請求項１ないし７のいずれかに記載の画像処理装置である。
【００２６】
この構成によれば、振幅画像データ生成手段によって生成される振幅画像データに基づいて表示装置における画像表示を行うことで、操作者に対して、個々の周波数成分の振幅分布を表す振幅画像を提供できる。この振幅画像は、たとえば、当該周波数成分の振幅データに対応する輝度を各画素に与えたものであってもよい。子宮蠕動に対応した周波数成分の振幅画像では、振幅の大きな画素が子宮の一部または全部の形状を形成する可能性が高い。したがって、このような振幅画像を提供することによって、操作者は、子宮蠕動に対応した周波数成分を容易に特定して選択することができる。
【００２７】
請求項９記載の発明は、前記フーリエ変換手段によって生成された少なくとも１つの周波数成分の位相データを用いて、子宮ＭＲ画像各部の当該周波数成分の位相分布を表す位相画像を表す位相画像データを生成する位相画像データ生成手段をさらに含むことを特徴とする請求項１ないし８のいずれかに記載の画像処理装置である。
この構成によれば、位相画像データ生成手段によって生成される位相画像データに基づいて表示装置における画像表示を行うことで、操作者に対して、個々の周波数成分の位相分布を表す位相画像を提供することができる。たとえば、この位相画像は、子宮ＭＲ画像の当該周波数成分の位相に対応した輝度を個々の画素に与えたものであってもよい。
【００２８】
子宮蠕動に対応した周波数成分について前述のような位相画像を作成すると、この位相画像には、子宮の一部または全部の形状に従って位相変化を示す画像部分が現れる。したがって、位相画像を提供することによって、子宮蠕動に対応する周波数成分を特定して選択する過程で、操作者を支援することができる。
或る周波数成分の位相画像を作成するときには、当該周波数成分の振幅が所定のしきい値に満たない画素については、画素値を固定（たとえば「０」に固定）することが好ましい。これにより、子宮蠕動に全く関係のないと考えられる部分の画像を潰すことができるので、操作者による判断がより容易になる。
【００２９】
請求項１０記載の発明は、前記画像データ記憶手段に記憶された子宮ＭＲ画像データ（複数枚の子宮ＭＲ画像データのなかから選択した代表的な１枚であってもよいし、複数の子宮ＭＲ画像データの平均値等をとった画像データであってもよい。）に前記位相勾配分布画像データ生成手段によって生成される位相勾配分布画像データをオーバレイして、オーバレイ画像データを生成する画像オーバレイ処理手段をさらに含むことを特徴とする請求項１ないし９のいずれかに記載の画像処理装置である。
【００３０】
この構成によれば、オーバレイ画像データに基づいて表示装置上に画像表示を行えば、子宮ＭＲ画像に位相勾配分布を表す着色画像がオーバレイ表示される。これにより、子宮の各部と、蠕動が生じている部位との対応関係を一目で把握できる画像を提供することができ、どの部位でどのような蠕動が生じているかを容易に判断できるようになる。
請求項１１記載の発明は、前記フーリエ変換手段によって生成された１つの周波数成分の位相データのうち、少なくとも２つの離間した画素の位相データに基づいて、当該２つの画素間の位相差を演算する位相差演算手段と、この位相差演算手段によって演算された位相差に基づいて、子宮の蠕動速度を計算する蠕動速度演算手段をさらに含むことを特徴とする請求項１ないし１０のいずれかに記載の画像処理装置である。
【００３１】
子宮蠕動に対応した周波数成分に関して２つの画素間の位相差を求めると、この位相差に基づいて子宮蠕動の速度を算出することができる。この算出された子宮蠕動の速度を、たとえば、表示装置に表示したりして提供することにより、子宮蠕動に関するより多くの情報を提供でき、より一層効果的に画像診断を支援できる。子宮蠕動の速度は、より具体的には、２つの画素間の位相差と、この２つの画素間の距離とに基づいて、算出することができる。２つの画素間の距離は、２つの画素の位置情報（座標）から算出できる。
【００３２】
請求項１２記載の発明は、操作者による画像中の２点の指示入力を受け付ける手段をさらに含み、前記位相差演算手段は、前記受け付けられた指示入力による２点の画素の位相データ（同じ周波数成分の位相データ）に基づいて、前記２点の画素間の位相差を演算するものであることを特徴とする請求項１１記載の画像処理装置である。
この構成によれば、入力装置から、速度の計測を行いたい２点を指示することができるので、操作者の所望する任意の２点間の蠕動速度の情報を提供することができる。
【００３３】
請求項１３記載の発明は、請求項１ないし１２のいずれかに記載の画像処理装置と、この画像処理装置に接続され、前記位相勾配分布画像を表示する表示装置とを含むことを特徴とする画像診断支援システムである。
この構成により、位相勾配分布画像データに対応した位相勾配分布画像を表示装置に表示して、この画像情報を使用者に提供することができる。これにより、子宮の蠕動の様子を明瞭に可視化した画像を提供でき、画像診断を良好に支援できる。
【００３４】
前記表示装置には、前述の最大振幅周波数分布画像、振幅画像および位相画像の表示も可能とされていることが好ましい。
請求項１４記載の発明は、操作者によって操作可能であり、前記画像処理装置に接続されて、この画像処理装置に対して操作者の指示入力を入力するための入力装置をさらに含むことを特徴とする請求項１３記載の画像診断支援システムである。
【００３５】
この構成により、使用者による指示入力によって、画像処理装置による処理内容を指示することができる。
請求項１５記載の発明は、子宮の蠕動運動を解析して可視化するための画像処理方法であって、子宮の蠕動の周期よりも短い周期で時系列に従って撮影した複数枚の子宮ＭＲ（磁気共鳴）画像に対応した複数枚の子宮ＭＲ画像データを準備するステップと、この複数枚の子宮ＭＲ画像データに対して時間軸方向のフーリエ変換を行い、子宮ＭＲ画像中の各画素について複数の周波数成分の位相データを生成するフーリエ変換ステップと、このフーリエ変換ステップによって生成された少なくとも１つの周波数成分の位相データを用いて、子宮ＭＲ画像各部における当該周波数成分の所定方向に関する位相勾配を算出し、その位相勾配を表す位相勾配データを生成する位相勾配データ生成ステップと、この位相勾配データ生成ステップによって生成された位相勾配データに基づき、画像各部を位相勾配に応じて着色した位相勾配分布画像を表す位相勾配分布画像データを、子宮の蠕動を表す画像データとして生成する位相勾配分布画像データ生成ステップとを含むことを特徴とする子宮蠕動可視化のための画像処理方法である。この方法により、請求項１の発明に関して説明した効果と同様の効果が得られる。
【００３６】
請求項１６記載の発明は、前記フーリエ変換ステップによって生成された１つの周波数成分の位相データのうち、少なくとも２つの離間した画素の位相データに基づいて、当該２つの画素間の位相差を演算する位相差演算ステップと、この位相差演算ステップによって演算された位相差に基づいて、子宮の蠕動速度を計算する蠕動速度演算ステップとをさらに含むことを特徴とする請求項１５記載の画像処理方法である。この方法により、請求項１１の発明と同様な効果が得られる。
【００３７】
むろん、画像処理方法の発明に関しても、画像処理装置の発明と同様な変形および改良が可能である。
請求項１７記載の発明は、子宮の蠕動の周期よりも短い周期で時系列に従って撮影した複数枚の子宮ＭＲ（磁気共鳴）画像に対応した複数枚の子宮ＭＲ画像データを処理し、子宮の蠕動運動を解析して可視化するための画像処理装置としてコンピュータを動作させるためのコンピュータプログラムであって、前記コンピュータを、前記複数枚の子宮ＭＲ画像データに対して時間軸方向のフーリエ変換を行い、子宮ＭＲ画像中の各画素について複数の周波数成分の位相データを生成するフーリエ変換手段、このフーリエ変換手段によって生成された少なくとも１つの周波数成分の位相データを用いて、子宮ＭＲ画像各部における当該周波数成分の所定方向に関する位相勾配を算出し、その位相勾配を表す位相勾配データを生成する位相勾配データ生成手段、およびこの位相勾配データ生成手段によって生成された位相勾配データに基づき、画像各部を位相勾配に応じて着色した位相勾配分布画像を表す位相勾配分布画像データを、子宮の蠕動を表す画像データとして生成する位相勾配分布画像データ生成手段として機能させることを特徴とする子宮蠕動可視化のためのコンピュータプログラムである。このコンピュータプログラムをコンピュータによって実行させることにより、このコンピュータを請求項１記載の画像処理装置として機能させることができる。
【００３８】
請求項１８記載の発明は、前記コンピュータを、さらに、前記フーリエ変換手段によって生成された１つの周波数成分の位相データのうち、少なくとも２つの離間した画素の位相データに基づいて、当該２つの画素間の位相差を演算する位相差演算手段、およびこの位相差演算手段によって演算された位相差に基づいて、子宮の蠕動速度を計算する蠕動速度演算手段として機能させることを特徴とする請求項１７記載のコンピュータプログラムである。このコンピュータプログラムをコンピュータに実行させることにより、このコンピュータを請求項１１の画像処理装置として機能させることができる。
【００３９】
むろん、コンピュータプログラムは、画像処理装置の発明に関連して説明した前述の種々の機能処理手段としてコンピュータを機能させるように変形または改良されてもよい。
請求項１９記載の発明は、前述のようなコンピュータプログラムを記録しており、コンピュータによる読み取りが可能な記録媒体である。このような記録媒体としては、ＣＤ−ＲＯＭおよびＤＶＤ−ＲＯＭに代表される光ディスク、ＭＯディスクに代表される磁気光ディスク、ならびにハードディスクドライブおよびフレキシブルディスクに代表される磁気ディスクを例示することができる。
【発明を実施するための最良の形態】
【００４０】
以下では、この発明の実施の形態を、添付図面を参照して詳細に説明する。
［子宮蠕動の解析の基本原理］
まず、この発明の一実施形態に係る画像処理装置による子宮ＭＲ画像の解析の基本原理について説明する。
図１に、子宮ＭＲ画像（矢状断面）の一例を示す。子宮は、領域ａに示す「体部」と、領域ｂに示す頚部とを有している。子宮の壁は、３層で構成されている。一番内側が内膜ｃであり、その内膜を覆うように存在するのが筋層ｄであり、内膜ｃと筋層ｄの境界がジャンクショナルゾーン（junctional zone）である。
【００４１】
図１に示す子宮ＭＲ画像は、高速スピンエコー法によって撮影された画像例である。スピンエコー法では、Ｔ１強調画像およびＴ２強調画像の撮影が可能であるが、Ｔ１強調画像では、子宮筋層と内膜が均一な低信号を持つため、これらを区別して撮影することができない。そのため、子宮のＭＲＩでは、Ｔ２強調画像での撮影が主として行われ、図１に示す子宮ＭＲ画像もＴ２強調画像である。
【００４２】
子宮の蠕動の解析には、子宮の蠕動の周期よりも短い時間間隔（たとえば２秒間隔）で時系列に従って撮影された複数枚の子宮ＭＲ画像のデータが用いられる。画像サイズは、たとえば、横１．２ｍｍ、縦１．４ｍｍである。撮影の性質上、１枚目の画像の輝度値は、他の画像に比べて高くなるので、１枚目は除き、２枚目からの画像を使用して解析を行うことが好ましい。後述する高速フーリエ変換のためには、２のべき乗枚（たとえば６４枚）の子宮ＭＲ画像が必要となる。したがって、たとえば、６０枚の子宮ＭＲ画像データが処理対象であり、２枚目以降の５９枚の子宮ＭＲ画像データに対して画像解析を行う場合には、これを、補間処理によって、６４枚の子宮ＭＲ画像データに変換して使用することになる。
【００４３】
複数枚の子宮ＭＲ画像をディスプレイ上に連続で表示して注意深く観察すると、子宮の内膜に内側に沿って微小な波状突起がほぼ一定の速さで移動し、子宮の内膜が波打つように運動をしていることが観察される。この微小波状突起の移動に伴い、内膜の画像部分を構成する各画素の輝度値が周期的な時間変化を示す。また、筋層では、黒い影のようなものが動いていることが観察できる。筋層は厚い筋肉から成り立っており、内膜と同様、周期的な運動をしていると考えられる。
【００４４】
図２に示す子宮ＭＲ画像において、子宮内膜上の点「１」、「２」および「３」の３点について、時間軸方向の輝度値の変化をプロットした例を図３に示す。横軸は、時系列に従う複数枚の画像の番号によって表される時間軸である。この図３から、子宮内膜の蠕動運動は、定点に注目したときに、濃度（輝度）変化の周期的な波の組合せとして表現できると推測できる。子宮筋層の蠕動運動についても同様である。
【００４５】
そこで、子宮ＭＲ画像の各画素の時系列データ（時系列に従って撮影された複数枚の子宮ＭＲ画像における同一位置の画素の輝度データ列）に対して、フーリエ変換処理が行われる。実際には、図２に示すように子宮部分以外の画像部分を排除するように矩形領域にトリミングした子宮ＭＲ画像を対象画像としてフーリエ変換処理が行われる。
ただし、子宮の画像部分の位置やその形状は、周囲の臓器の運動等の影響を受けて変動する。そのため、処理領域を固定して解析を行うフーリエ変換では、子宮の位置が変化すると、運動解析の対象位置がずれてしまい、適切な解析を行うことができない。そこで、フーリエ変換に先だって、複数枚の子宮ＭＲ画像間の子宮の位置および形状の変動を補償するための幾何学変換処理が行われる。具体的には、画像間の変形を近似するアフィン変換パラメータを非線形最適化を用いて求め、そのパラメータを用いたアフィン変換によって、子宮の位置および形状が合わせ込まれる。より詳細には、まず並進移動成分を決定した後、その並進移動のパラメータを初期値としてアフィン変換パラメータが求められる。このアフィン変換パラメータを用いて原画像が変形され、対応する領域の切り出し（トリミング）を行った画像列を作成し、その画像列に対してフーリエ変換が行われる。
【００４６】
フーリエ変換には、離散フーリエ変換の一形態である高速フーリエ変換（ＦＦＴ：Fast Fourier Transform）が適用される。これにより、各画素の輝度変化における主な周波数成分の振幅および位相が求められる。これに基づいて、子宮内膜および筋層の蠕動の方向および速度を解析することができる。
より具体的に説明すると、Ｎ枚の子宮ＭＲ画像中の同じ位置にある１つの画素についてのＮ個の時系列データ（輝度値を表す画像データ）をｘ0，ｘ1，ｘ2，…，ｘn，…，ｘN-1と表すと、周波数ωk＝２πｋ／Ｎ（ｋ＝０，１，２，３，……，Ｎ−１）についての離散フーリエ変換Ｘkは、次の式(1)で与えられる。ただし、ｊは虚数単位である。
【００４７】
【数１】
【００４８】
そして、周波数ωkの周波数成分の振幅データ｜Ｃk｜および位相データφkは、それぞれ、次の式(2)および式(3)で与えられる。
｜Ｃk｜＝｛ＸkRe2＋ＸkIm2｝1/2 ……（２）
φk ＝ｔａｎ-1（−ＸkIm／ＸkRe） ……（３）
ただし、ＸkReおよびＸkImは、それぞれ、Ｘkの実数部および虚数部である。
【００４９】
これらの振幅データ｜Ｃk｜および位相データφkが、子宮ＭＲ画像の個々の画素について、計算可能なすべての周波数ωkの周波数成分に関して求められる。
高速フーリエ変換を行うためには、前述のように、２のべき乗個の時系列データが必要となる。したがって、有効な子宮ＭＲ画像が５９枚であるときは、この５９枚の子宮ＭＲ画像（アフィン変換後に切り出した画像）の個々の画素の輝度値を５９個から６４個に補間したデータを予め作成し、この６４個のデータに対して高速フーリエ変換が行われる。
【００５０】
この高速フーリエ変換によって得られた振幅データに基づき、周波数「０」（直流成分）以外で最大振幅をとる周波数が、各画素について求められる。そして、各画素の最大振幅の周波数成分を表す最大振幅周波数分布画像がディプレイ上に表示される。
具体的には、図４(a)に示すように、複数の周波数成分に対して予め異なる色が割り当てられ、個々の画素が、当該画素の最大振幅の周波数成分に応じて着色される。こうして、最大振幅の周波数成分分布が、ディプレイ上に疑似カラー表示されることになる。
【００５１】
子宮の蠕動運動は低周波成分に現れるため、図４(a)の例では、第１周波数〜第１０周波数（ｋ＝１〜１０）の周波数成分についての分布画像を求めた例が示されている。第１周波数とは、６４枚で１周期となる周波数成分であり、第２周波数とは、３２（＝６４／２）枚で１周期となる周波数成分である。以下、第３周波数成分以上の周波数成分についても同様である。
【００５２】
参考のため、図４(b)には図４(a)のカラー画像中の第１周波数の色部分を強調した画像を示し、図４(c)には図４(a)のカラー画像中の第２周波数の色部分を強調した画像を示し、図４(d)には図４(a)のカラー画像中の第６周波数の色部分を強調した画像を示す。
内膜部分または筋層部分の蠕動は、規則的な周期運動であるので、内膜部分または筋層部分の画素値（最大振幅の周波数成分）は同じであることが期待される。実際、図４の例では、第６周波数成分を最大振幅周波数成分とする複数の画素が、子宮内膜に近似した形状の画像部分を形成している。
【００５３】
したがって、最大振幅周波数分布画像を見ることによって、使用者（操作者。一般的には専門医）は、子宮の蠕動に対応した周波数成分を容易に特定することができる。
子宮蠕動に対応した周波数成分の特定を補助する目的で、必要に応じて、振幅画像データおよび位相画像データが作成され、それらに対応する振幅画像および位相画像がそれぞれディスプレイ上に表示される。
【００５４】
振幅画像データは、一つの周波数成分について、各画素の振幅データを画素値（輝度値）とした振幅分布データである。このような振幅画像データに対応する振幅画像は、当該周波数成分の振幅が大きい画素ほど明るく表示される。子宮の蠕動は、前述のとおり、規則的な周期運動であるから、子宮の蠕動に対応した周波数成分についての振幅画像は、子宮の内膜部分または筋層部分に近似した形状の画像部分を形成すると期待される。したがって、振幅画像中に子宮の形状に近似した画像部分が形成される周波数成分を、子宮蠕動に対応するものと特定することができる。振幅画像の一例は、図１６に示されている。
【００５５】
位相画像データは、一つの周波数成分について、各画素の位相データを画素値（輝度）とした位相分布データである。ただし、当該周波数成分についての振幅データが所定のしきい値未満の画素については、画素値を「０」に固定する。このような位相画像データは、位相が大きい画素ほど明るく表示される。しきい値未満の画素の画素値を「０」とするので、子宮蠕動に対応した周波数成分の位相画像中には、やはり、子宮の形状に対応した画像部分が形成されることが期待でき、さらに、この画像部分は、子宮の蠕動方向および速度に対応した輝度変化を示すと期待される。したがって、このような特徴の画像部分が位相画像中に見られる周波数成分を、子宮蠕動に対応するものと特定することができる。位相画像の一例は、図１７に示されている。
【００５６】
こうして子宮蠕動に対応した周波数成分が特定されると、位相画像の濃度勾配（輝度勾配）が算出される。すなわち、個々の画素について、所定方向に関して、隣接画素との画素値の差分（濃度勾配）が求められる。この差分は、結局、位相勾配を表す位相勾配データにほかならない。位相勾配を求める方向は、たとえば、操作者がディスプレイに表示された位相画像上で指定する。
【００５７】
さらに、位相勾配データの値に応じて色情報が割り当てられる。たとえば、位相勾配データが正であれば、第１の色（たとえば青色）が当該画素に割り当てられ、位相勾配データが負であれば、第１の色とは異なる第２の色（たとえば赤色）が当該画素に割り当てられる。これにより、子宮蠕動の方向を色情報で表した位相勾配分布画像が形成される。位相勾配分布画像の例は、図２２、図２８および図３６に示されている。
【００５８】
むろん、位相勾配の値に応じて画素値（輝度）を異ならせてもよい。これにより、子宮蠕動の速さを併せて表示することができる。
子宮蠕動の速度は、画像中の任意の２点を指定することによって行える。速度＝距離÷時間であり、時間＝周期×位相差÷２πである。「周期」は１周期の時間であり、周期＝画像枚数×２（秒）÷周波数で表される。周波数とは、子宮蠕動に対応すると特定された周波数成分の周波数である。
【００５９】
したがって、子宮蠕動に対応すると特定された周波数成分について、指定された２点の位相データを用いて位相差を演算すれば、この位相差に基づいて、子宮蠕動の速度を求めることができる。
［画像診断支援システムの構成］
図５は、この発明の一実施形態に係る画像診断支援システムの基本的なハードウェア構成を説明するためのブロック図である。この画像診断支援システムは、コンピュータシステムとしての基本構成を有している。すなわち、この画像診断支援システムは、コンピュータ本体１０と、このコンピュータ本体１０に接続された入力装置１１としてのキーボード１１Ａおよびマウス（ポインティングデバイスの一例）１１Ｂと、コンピュータ本体１０に接続されたディプレイ（表示装置）１３と、コンピュータ本体１０に接続されたプリンタ１４とを備えている。
【００６０】
コンピュータ本体１０は、演算処理部としてのＣＰＵ１５と、主記憶装置としてのＲＯＭ１６およびＲＡＭ１７と、補助記憶装置してのハードディスクドライブ（ＨＤＤ）１８およびＣＤ−ＲＯＭドライブ１９と、入力制御部２０と、ディスプレイインタフェース２１と、プリンタインタフェース２２とを備えており、これらはバス２３によって互いに接続されている。入力制御部２０は、入力装置１１からの入力信号を処理してＣＰＵ１５に入力する。ディスプレイインタフェース２１は、画像データの入力を受けて、その画像データに対応する画像をディスプレイ１３に表示させる働きを担う。プリンタインタフェース２２は、プリンタ１４との接続のためのものである。
【００６１】
ＣＤ−ＲＯＭドライブ１９には、コンピュータプログラムを記録した記録媒体としてのＣＤ−ＲＯＭ（ディスク）２５を装填することができる。画像処理プログラムを記録したＣＤ−ＲＯＭ２５をＣＤ−ＲＯＭドライブ１９に装填し、必要なインストール操作を行ったうえで、当該画像処理プログラムをコンピュータ本体１０によって実行させることにより、このコンピュータ本体１０は画像処理装置として機能し、当該コンピュータシステムは画像診断支援システムとして機能することができる。
【００６２】
むろん、コンピュータプログラムは、記録媒体に記録された形態で提供される必要はなく、ネットワークを介して提供されてもよい。すなわち、コンピュータ本体１０にネットワークインタフェースを備え、このネットワークインタフェースを介してコンピュータプログラムをコンピュータ本体１０にダウンロードすることによっても、コンピュータ本体１０を画像処理装置として機能させることができる。
【００６３】
一方、処理対象の子宮ＭＲ画像データは、たとえば、ハードディスクドライブ１８に格納される。この子宮ＭＲ画像データも、ＣＤ−ＲＯＭドライブ１９を経由してハードディスクドライブ１８に格納されてもよいし、ネットワークを介して取得してハードディスクドライブ１８に格納するようにしてもよい。子宮ＭＲ画像データは、子宮の蠕動運動の周期よりも短い時間間隔（たとえば２秒）で時系列に従って撮影された複数枚の子宮ＭＲ画像を表すデータである。
【００６４】
図６は、前記画像処理装置として機能するコンピュータ本体１０の機能的な構成を説明するためのブロック図である。コンピュータ本体１０は、ＣＰＵ１５によって画像処理プログラムを実行させることによって、図６に示す機能処理部を実質的に有することになる。すなわち、コンピュータ本体１０は、画像処理演算部３０と、データ記憶部５０と、画像表示制御部７０とを有している。画像処理演算部３０および画像表示制御部７０は、主として、ＣＰＵ１５による演算処理によって実現される機能処理部である。また、データ記憶部５０は、ＲＡＭ１７およびハードディスクドライブ１８の記憶領域で構成される。
【００６５】
画像処理演算部３０は、アフィン変換処理部３１と、トリミング処理部３２と、フーリエ変換処理部３３と、最大振幅周波数成分抽出処理部３４と、最大振幅周波数分布画像データ生成処理部３５と、周波数成分選択入力受付け部３６と、振幅画像データ生成処理部３７と、位相画像データ生成処理部３８と、位相勾配算出方向入力受付け部３９と、位相勾配データ生成処理部４０と、位相勾配分布画像データ生成処理部４１と、画像オーバレイ処理部４２と、速度測定点入力受付け部４３と、位相差演算部４４と、蠕動速度演算部４５とを備えている。
【００６６】
また、データ記憶部５０は、原子宮ＭＲ画像データ記憶部５１と、アフィン変換後子宮ＭＲ画像データ記憶部５２と、対象子宮ＭＲ画像データ記憶部５３と、フーリエ変換データ記憶部５４と、最大振幅周波数成分データ記憶部５５と、最大振幅周波数分布画像データ記憶部５６と、振幅画像データ記憶部５７と、位相画像データ記憶部５８と、位相勾配データ記憶部５９と、位相勾配分布画像データ記憶部６０と、オーバレイ画像データ記憶部６１と、位相差データ記憶部６２と、蠕動速度データ記憶部６３とを備えている。
【００６７】
原子宮ＭＲ画像データ記憶部５１は、ＭＲＩによって子宮の蠕動運動の周期よりも短い時間間隔（たとえば２秒）で時系列に従って撮影された複数枚の子宮ＭＲ画像（原データ。ただし、必要に応じて補間処理をしたもの）を記憶する。
アフィン変換処理部３１は、原子宮ＭＲ画像データ記憶部５１から子宮ＭＲ画像データを読み出し、この子宮ＭＲ画像データにアフィン変換処理を行い、アフィン変換後の子宮ＭＲ画像データをアフィン変換後子宮ＭＲ画像データ記憶部５２に格納する。
【００６８】
トリミング処理部３２は、アフィン変換処理後の子宮ＭＲ画像データをアフィン変換後子宮ＭＲ画像データ記憶部５２から読み出し、この子宮ＭＲ画像データから子宮部分以外の画像を排除するトリミング処理を行い、この処理後のトリミングされた子宮ＭＲ画像データを対象子宮ＭＲ画像データ記憶部５３に格納する。
フーリエ変換処理部３３は、このトリミング処理部３２による処理後の子宮ＭＲ画像データを対象子宮ＭＲ画像データ記憶部５３から読み出し、この子宮ＭＲ画像データに対してフーリエ変換（高速フーリエ変換）を行う。より具体的には、フーリエ変換処理部３３は、アフィン変換を経てトリミングされた複数枚の子宮ＭＲ画像データに対して、時間軸方向のフーリエ変換を行う。これにより、子宮ＭＲ画像の個々の画素について、複数の周波数成分に関して、振幅データおよび位相データが求められる。すなわち、子宮ＭＲ画像の各画素の輝度値の時間変化が複数の周波数成分に分解され、各周波数成分について、その振幅を表す振幅データおよびその位相を表す位相データが求められる。これらのデータは、フーリエ変換データ記憶部５４に格納される。
【００６９】
最大振幅周波数成分抽出処理部３４は、フーリエ変換データ記憶部５４に格納された各画素の複数の周波数成分についての振幅データに基づき、子宮ＭＲ画像の個々の画素について、最大の振幅データを有する周波数成分を抽出する。こうして、画素ごとに抽出された周波数成分を表すデータは、最大振幅周波数成分データ記憶部５５に記憶される。
最大振幅周波数分布画像データ生成処理部３５は、最大振幅周波数成分データ記憶部５５に記憶されたデータに基づき、前記抽出された最大振幅周波数成分の分布を表す画像（最大振幅周波数分布画像）を表す最大振幅周波数分布画像データを生成し、最大振幅周波数分布画像データ記憶部５６に格納する。具体的には、最大振幅周波数分布画像データ生成処理部３５は、複数の周波数成分に対して異なる色データを割り当て、当該色データを画素値とする最大振幅周波数分布画像データを生成する。画像表示制御部７０は、この最大振幅周波数分布画像データを読み出して、ディスプレイ１３（図５参照）に最大振幅周波数分布画像を表示させる。これにより、子宮ＭＲ画像中における最大振幅周波数成分の分布を表す疑似カラー表示画像がディスプレイ１３上に表示される。
【００７０】
周波数成分選択入力受付け部３６は、入力装置１１からの周波数成分選択入力を受け付ける。操作者は、ディスプレイ１３に表示された最大振幅周波数分布画像から判断して、子宮蠕動に対応した周波数成分を特定し、この周波数成分の選択指示を入力装置１１を介して入力する。この選択指示が、周波数成分選択入力受付け部３６によって受け付けられる。
【００７１】
振幅画像データ生成処理部３７は、周波数成分選択入力受付け部３６によって受け付けられた選択指示に対応する周波数成分についての振幅画像の表示のための振幅画像データを生成する。具体的には、振幅画像データ生成処理部３７は、フーリエ変換データ記憶部５４から、前記選択された周波数成分についての各画素の振幅データを読み出し、この振幅データを画素値とする振幅画像データ（振幅分布データ）を作成し、この振幅画像データを振幅画像データ記憶部５７に格納する。画像表示制御部７０は、この振幅画像データを読み出して、ディスプレイ１３の表示制御を行う。これにより、前記選択された周波数成分についての振幅データを各画素の輝度値とする振幅画像がディスプレイ１３上に表示される。
【００７２】
一方、位相画像データ生成処理部３８は、周波数成分選択入力受付け部３６によって受け付けられた選択指示に対応する周波数成分についての位相画像の表示のための位相画像データを生成する。具体的には、位相画像データ生成処理部３８は、フーリエ変換データ記憶部５４から、前記選択された周波数成分についての各画素の位相データを読み出し、この位相データを画素値とする位相画像データ（位相分布データ）を作成し、この位相画像データを位相画像データ記憶部５８に格納する。ただし、位相画像データ生成処理部３８は、振幅画像データ記憶部５７に記憶された振幅画像データを参照し、その値（すなわち振幅データの値）が、所定のしきい値に達していない画素については、位相画像データを「０」に固定する。画像表示制御部７０は、位相画像データ記憶部５８から位相画像データを読み出して、ディスプレイ１３の表示制御を行う。これにより、位相画像データが「０」に固定された画素以外の画素について、前記選択された周波数成分についての位相データを各画素の輝度値とする位相画像がディスプレイ１３上に表示されることになる。
【００７３】
位相勾配算出方向入力受付け部３９は、入力装置１１からの位相勾配算出方向指示入力を受け付ける。操作者は、ディスプレイ１３に表示される振幅画像および／または位相画像から、子宮蠕動の方向を判断し、その方向を入力装置１１から指示する。この指示入力が位相勾配算出方向入力受付け部３９によって受け付けられることになる。
位相勾配データ生成処理部４０は、位相画像データ記憶部５８に記憶された位相画像データを読み出し、個々の画素について、前記位相勾配算出方向入力受付け部３９によって受け付けられた位相勾配算出方向に隣接する画素との間での差分値（微分値）を求め、この差分値を位相勾配データとして生成し、位相勾配データ記憶部５９に格納する。位相勾配データは、したがって、子宮ＭＲ画像の各部における運動の速度（方向および速さ）を表す。位相画像データ記憶部５８に、子宮蠕動に対応した周波数成分の位相画像データが格納されていると、位相勾配データは、子宮の各部の蠕動の速度（方向および速さ）を表すことになる。
【００７４】
位相勾配分布画像データ生成処理部４１は、位相勾配データ記憶部５９から位相勾配データを読み出し、位相勾配分布をカラー表示するための位相勾配分布画像データを生成する。具体的には、位相勾配分布画像データ生成処理部４１は、位相勾配データの値ごとに予め異なる色データを割り当てており、このような色データを画素値とする画像データを位相勾配分布画像データとして生成し、位相勾配分布画像データ記憶部６０に格納する。画像表示制御部７０は、位相勾配分布画像データ記憶部６０から位相勾配分布画像データを読み出し、ディスプレイ１３の表示制御処理を実行する。これにより、位相勾配分布を疑似カラー表示した位相勾配分布画像がディスプレイ１３上に表示されることになる。
【００７５】
たとえば、位相勾配分布画像データ生成処理部４１は、位相勾配データの値が正であれば第１の色データ（たとえば青色データ）を位相勾配分布画像データとし、位相勾配データの値が負であれば第１の色データとは異なる第２の色データ（たとえば赤色データ）を位相勾配分布画像データとするものであってもよい。この場合、入力装置１１から、子宮頚部から子宮体部に向かう方向（またはその正反対方向）を位相勾配算出方向として指示しておけば、頚部から体部に向かう蠕動と、体部から頚部に向かう蠕動とを、色分けして表示することができる。
【００７６】
また、位相勾配分布画像データ生成処理部４１は、位相勾配データの値が正であれば第１の色データ（たとえば青色データ）を割り当て、位相勾配データの値が負であれば第１の色データとは異なる第２の色データ（たとえば赤色データ）を割り当てるとともに、さらに、位相勾配データの絶対値を輝度値（画素値）とする位相勾配分布画像データを生成するものであってもよい。このようにすれば、各部の蠕動の方向を色分けして表すとともに、各部の蠕動の速さを輝度によって表すことができる。
【００７７】
画像オーバレイ処理部４２は、原子宮ＭＲ画像データ記憶部５１から子宮ＭＲ画像データ（たとえば、アフィン変換処理の際の基準とした１枚の画像）を読み出し、位相勾配分布画像データ記憶部６０から位相勾配分布画像データを読み出し、これらをオーバレイしたオーバレイ画像データを生成して、オーバレイ画像データ記憶部６１に格納する。すなわち、画像オーバレイ処理部４２は、子宮ＭＲ画像データにおいて、トリミング処理部３２によって切り出された矩形領域に対応する領域に、位相勾配分布画像データを重ね合わせる。画像表示制御部７０は、オーバレイ画像データ記憶部６１からオーバレイ画像データを読み出し、これをディスプレイ１３に画像表示する。これによって、子宮ＭＲ画像（原画像）に位相勾配分布画像が重ね合わせられて表示されるので、子宮の蠕動の様子を、その位置とともに一目で把握できるようになる。
【００７８】
むろん、画像表示制御部７０により、位相勾配分布画像データ記憶部６０に記憶された画像データをディスプレイ１３に画像表示させることもできる。このときには、位相勾配分布画像のみがディスプレイ１３に表示される。
速度測定点入力受付け部４３は、位相画像、位相勾配分布画像またはオーバレイ画像がディスプレイ１３に表示されている状態で、入力装置１１からの画像中の２点の入力を受け付ける。この入力は、当該２点間の速度を計測すべき速度測定点として受け付けられる。
【００７９】
位相差演算部４４は、速度測定点入力受付け部４３によって受け付けられた２点の画素の位相データを位相画像データ記憶部５８から読み出し、それらの差を位相差データとして演算して、位相差データ記憶部６２に格納する。
蠕動速度演算部４５は、速度測定点入力受付け部４３によって受け付けられた２つの速度測定点間の距離と、位相差データ記憶部６２に記憶された位相差データとに基づき、前記２点間の速度を蠕動速度として求め、この蠕動速度を表す蠕動速度データを蠕動速度データ記憶部６３に格納する。
【００８０】
画像表示制御部７０は、蠕動速度データ記憶部６３に格納された蠕動速度データをディスプレイ１３上に表示（たとえば数値表示）する。
図７は、子宮ＭＲ画像データを処理して子宮の蠕動を可視化するための処理の流れをまとめて示すフローチャートである。ＭＲＩによって得られた複数枚の子宮ＭＲ画像データに対して、画像間の子宮の位置および形状の変動を補償するための幾何学変換であるアフィン変換処理が行われ（ステップＳ１）、さらに、子宮部分の画像を切り出すためのトリミングが行われる（ステップＳ２）。次に、トリミングされた複数枚の子宮ＭＲ画像データを対象画像データとして以下の処理が行われる。
【００８１】
すなわち、複数枚の子宮ＭＲ画像データに対して、時間軸方向のフーリエ変換処理（ＦＦＴ）が実行され、個々の画素について、複数の周波数成分に関する振幅データおよび位相データが生成される（ステップＳ３）。次いで、個々の画素に関して、最大の振幅データを有する周波数成分（最大振幅周波数成分）が抽出され（ステップＳ４）、この最大振幅周波数成分に対応した色データを画素値とした疑似カラー画像である最大振幅周波数分布画像データが作成され、これに対応する画像（最大振幅周波数分布画像）がディスプレイ１３に表示される（ステップＳ５）。
【００８２】
操作者（典型的には、画像診断を行う専門医）は、ディスプレイ１３に表示された最大振幅周波数分布画像に基づき、子宮の蠕動に対応すると考えられる周波数成分を選択して、選択指示を入力装置１１から入力する（ステップＳ６）。
これにより、選択された周波数成分の振幅画像データおよび位相画像データが作成され、これらにそれぞれ対応する振幅画像および位相画像がディスプレイ１３に表示される（ステップＳ７，Ｓ８）。操作者は、表示された振幅画像および位相画像から、選択した周波数成分が子宮蠕動に対応しているかどうかを判断する（ステップＳ９）。当該選択した周波数が子宮蠕動を反映していると考えられる場合には、操作者は、以降の処理の実行を、入力装置１１から指示する（ステップＳ１０）。これにより、ステップＳ１１以降の処理が実行される。一方、前記選択された周波数成分が、子宮蠕動を反映していないと考えられる場合には、操作者は、別の周波数成分の選択指示を入力装置１１から入力する（ステップＳ６）。これにより、新たに選択された当該別の周波数成分に関して、振幅画像データおよび位相画像データが作成され、対応する振幅画像および位相画像がディスプレイ１３にそれぞれ表示される（ステップＳ７，Ｓ８）。
【００８３】
ステップＳ１１では、操作者は、子宮蠕動を反映していると判断された周波数成分の位相画像上の２点を入力装置１１から指定することにより、位相勾配の算出方向を指定する。これにより、当該方向に沿って、位相画像中の画素間の輝度勾配が位相勾配データとして求められる（ステップＳ１２）。さらに、この位相勾配データの符号に対応した２色の色データを画素値とする位相勾配分布画像データが作成され、対応する位相勾配分布画像がディスプレイ１３に表示される（ステップＳ１３）。さらに、必要に応じて、使用者は、子宮ＭＲ画像データと位相画像分布画像とのオーバレイ画像の表示を、入力装置１１から指示する（ステップＳ１４）。これに応答して、両画像をオーバレイしたオーバレイ画像データが作成されて、対応するオーバレイ画像がディスプレイ１３上に表示される（ステップＳ１５）。
【００８４】
子宮蠕動の速度を求める場合には、位相画像（ステップＳ８）、位相勾配分布画像（ステップＳ１３）およびオーバレイ画像（ステップＳ１５）のいずれかにおいて、操作者は、入力装置１１から画像中の２点を速度測定点として入力する（ステップＳ１６）。これに応答して、入力された２つの速度測定点間の位相データの差（位相差）が求められ（ステップＳ１７）、この位相差を用いて、蠕動運動の速度が求められ（ステップＳ１８）、そして、求められた蠕動運動の速度がディスプレイ１３上に表示される（ステップＳ１９）。
［シミュレーション実験］
子宮ＭＲ画像に人工的に内膜運動を付け加えた画像を作成し、その画像に対して前述の画像処理を行った。
【００８５】
図８は、シミュレーション実験のために作成した６４枚の時系列子宮ＭＲ画像を説明するための画像例である。１枚の対象画像（子宮ＭＲ画像）上の位置Ａに黒い点を付け加えて１枚目の子宮ＭＲ画像とし、２枚目以降の画像には同じ対象画像上に、順次、位置Ｂに向かって位置をずらせながら同じく黒い点を付け加え、８枚目で黒い点が位置Ｂに到達するようにした。９枚目以降の画像についても同様に、同じ対象画像を用い、８枚で位置Ａから位置Ｂまで移動するように黒い点を画像中に付け加えた。こうして、８枚で１周期となる模擬蠕動運動を付け加えた６４枚の子宮ＭＲ画像データを作成した。
【００８６】
位置ＡおよびＢの画素の輝度値の時間変化を図９に示す。横軸は、時間であるが、１枚目から６４枚目までの画像の番号で表されている。この図より、黒い点が移動することによって、位置Ａ，Ｂの画素の輝度値が変化している様子が観測できる。黒い点が位置Ａ，Ｂに位置する時間が異なるので、位置Ａ，Ｂの画素にそれぞれ対応した２つの曲線（折れ線）が、時間軸方向にずれている。
【００８７】
図８の例において、位置Ａの座標は（１１９，８７）、位置Ｂの座標は（９１，１１１）であった。この場合、２点間の距離は４４．６５（ｍｍ）となる。つまり、（８−１）×２（ｓｅｃ）で４４．６５（ｍｍ）を移動するため、ここでシュミレーションデータとして作成した運動の速度は、４４．６５（ｍｍ）／１４（ｓｅｃ）＝３．１９（ｍｍ／ｓｅｃ）となる。
【００８８】
前述のような６４枚の子宮ＭＲ画像データに対して時間軸方向にフーリエ変換を施し、各画素の最大振幅の周波数成分を求め、最大振幅周波数分布画像を作成した。これを図１０（(a)はカラー画像、(b)は(a)のカラー画像の第８周波数の色部分を強調した画像）に示す。この最大周波数分布画像から、第８周波数成分（＝６４／８）が運動を表現していることが分かる。この第８周波数成分に対する振幅画像および位相画像を図１１および図１２にそれぞれ示す。
【００８９】
振幅画像（図１１）において、高輝度となっている場所が、人工画像（黒い点を付加した６４枚の子宮ＭＲ画像）で黒い点が移動した場所である。位相画像（図１２）では、所定のしきい値未満の振幅の画素の画素値は「０」に固定されるので、図１１の振幅画像において高輝度となっている領域と同じ領域に有意な画素値を持つ画素が現れている。この位相画像の画素値（位相の値）を用いることで、２点間の運動の速度（向きと速さ）を解析できる。つまり、この場合、高輝度な位置Ａの方から低輝度な位置Ｂに向かう運動として検出することができる。
【００９０】
位相画像中の位置Ａと位置Ｂは、人工画像を作成した際の運動の始点と終点である。この２点の距離と位相差より、運動の速度を求める。位置（１１９，８７）と位置Ｂ（９１，１０８）の間の距離は４４．６５（ｍｍ）であり、位相はそれぞれ３．１４（ｒａｄ），−２．３６（ｒａｄ）であった。したがって、位置Ａ，Ｂの画素間の位相差は５．５０（ｒａｄ）である。よって、運動の速度は、３．１９（ｍｍ／ｓｅｃ）となり、設定した速度に等しい結果が得られた。
【００９１】
位相勾配を求める方向を、位置Ａから位置Ｂに向かう方向に指定して、各画素ごとに隣接画素との輝度値の差（位相勾配）を求め、正の位相勾配（位置Ａから位置Ｂに向かって位相が増加する場合に対応）に対して青色の色データを対応付け、負の位相勾配（位置Ｂから位置Ａに向かって位相が増加する場合に対応）に対して赤色の色データを対応付け、これらの色データを画素値とする位相勾配分布画像を作成すると、図１３に示す画像が得られる。すなわち、位置Ａ，Ｂ間の領域が青色に着色された画像が得られる。
［実験１（排卵期の画像）］
図１４に示す対象画像（６４枚の子宮ＭＲ画像。図１４にはその中の代表的な１枚を示す。）は或る被験者の排卵期の子宮をＭＲＩによって撮影したものである。この対象画像に本実施形態の画像処理を適用した。
【００９２】
図１５(a)に最大振幅周波数分布画像（カラー画像）を示す。図１５(b)には図１５(a)のカラー画像中の第３周波数の色部分を強調した画像を示し、図１５(c)には図１５(a)のカラー画像中の第６周波数の色部分を強調した画像を示し、図１５(d)には図１５(a)のカラー画像中の第７周波数の色部分を強調した画像を示す。
図１５(a)の画像から、子宮内膜の運動の周波数が一定でないこと、すなわち、内膜の運動にはいくつかの異なる周波数の運動があることが分かる。この例では、最大振幅周波数分布画像から一意に適切な周波数成分（内膜の蠕動に対応する周波数成分）を特定できない。そこで、振幅画像と位相画像をいくつかの周波数成分について作成した。その結果、第７周波数成分の振幅画像および位相画像がそれぞれ図１６および図１７に示すとおりとなり、子宮内膜の形状の画像部分が形成された。したがって、第７周波数成分が子宮蠕動に対応するものであると判断した。
【００９３】
図１７において、内膜部分が、子宮の頚部から体部に向かって徐々に暗くなっているため、対象画像子宮の内膜の運動方向は、頚部から体部に向かう方向であることが分かる。
図１７の４つの点ａ，ｂ，ｃ，ｄの輝度値の時間軸方向の変化を図１８および図１９に示す。図１８は原データ（フーリエ変換前のデータ）、図１９は直流成分と第７周波数成分に逆フーリエ変換を適用したものである。図１９を見ると、４つの点ａ，ｂ，ｃ，ｄの輝度値変化が時間軸方向にずれていることが分かる。この時間軸方向のずれが、図１７の位相画像に輝度値の変化として現れている。
【００９４】
排卵期には、精子の輸送のため、子宮の内膜が頚部から体部へ運動していると考えられており、排卵期の子宮の画像を扱った本実験では、そのことが確かめられた。
内膜上部の速度は図１７中の点ａと点ｄの位相から求めた。点ａの座標は（１１４，７８）で、点ｄの座標は（７５，９１）であったので、２点間の距離は５０．２１（ｍｍ）となる。位相はそれぞれ１．８９（ｒａｄ）、−１．５７（ｒａｄ）であったので、位相差は３．４６（ｒａｄ）である、よって、速度は５．３２（ｍｍ／ｓｅｃ）となる。内膜下部の速度は、図１７の画像中の点６と点７から求めた。点６の座標は（１０９，９４）で、点７の座標は（８９，１０９）であったので、２点間の距離は３１．８９（ｍｍ）となる。位相はそれぞれ１．０９（ｒａｄ），−１．７１（ｒａｄ）だったので、位相差は２．８０（ｒａｄ）である。よって、速度は４．１７（ｍｍ／ｓｅｃ）となる。
【００９５】
本手法により解析された運動が、画像列を連続で表示したときに視覚で捉えられる運動と同じ運動を捉えているかどうかを評価するために、画像中で観察できる運動の速度および位相から運動の速度を求めた。
図２０は、前記の実験の対象画像の画像列の２８枚目の子宮ＭＲ画像を示し、図２１は当該画像列の３７枚目の子宮ＭＲ画像を示す。この実験の対象画像では、図２０に示す２８枚目から図２１に示す３７枚目の画像の間で、黒い影（図中の矢印で示す部分）が、内
膜上部では点Ｃ（１２３，７３）から点Ｄ（５７，１０３）まで８９．６５（ｍｍ）移動し、内膜下部では点Ｅ（１１９，８７）から点Ｆ（７１，１１５）まで６９．６７ｍｍ移動していた。経過時間は１８（ｓｅｃ）（＝２sec×９）であるため、内膜上部と下部の運動の速度は、それぞれ４．９８（ｍｍ／ｓｅｃ）、３．８７（ｍｍ／ｓｅｃ）であった。 下記表１に示すように、本手法により求めた速度と、画像中で捉えられる運動の速度がそれぞれほぼ近い値となっており、本手法で解析した結果と、目視で検出した結果とが、同等な運動を捉えていると考えられる。
【００９６】
【表１】
【００９７】
よって、本手法が、一般には手作業で計測することが困難な子宮の蠕動運動を、画像解析によって求めることができる有効な手法であることが確認された。
位相勾配を求める方向を、図１７の点ａから点ｄ向かう方向に指定して、各画素ごとに隣接画素との輝度値の差（位相勾配）を求め、正の位相勾配（点ａから点ｄに向かって位相が増加する場合に対応）に対して青色の色データを対応付け、負の位相勾配（点ｄから点ａに向かって位相が増加する場合に対応）に対して赤色の色データを対応付け、これらの色データを画素値とする位相勾配分布画像を作成すると、図２２（(a)はカラー画像、(b)は(a)のカラー画像の青色部分を強調した図。右下部の矢印は、左下に向かうものが青色矢印であり、右上に向かうものが赤色矢印である。）に示す画像が得られる。すなわち、子宮内膜に対応する部分が青色に着色された画像が得られ、子宮頚部から体部に向かう蠕動運動が明瞭に可視化されている。
【００９８】
このカラー画像を原データに重ね合わせることにより、図２３に示すオーバレイ画像が得られる（(a)はカラー画像、(b)は(a)のカラー画像の青色領域を白抜きで強調表示した画像）。
［実験２（月経期の画像］
図２４に示す対象画像（６４枚の子宮ＭＲ画像。図２４にはそのなかの代表的な１枚を示す。）は、実験１の場合と同じ被験者の月経期の画像である。この画像に前記実施形態の手法を適用し、排卵期と月経期の運動の違いを検証した。
【００９９】
図２５(a)に最大振幅周波数分布画像を示す。参考のため、図２５(b)には図２５(a)のカラー画像中の第２周波数の色部分を強調した画像を示し、図２５(c)には図２５(a)のカラー画像中の第３周波数の色部分を強調した画像を示す。
図２５(a)の画像から、主な運動の周波数は第２周波数成分であると判断した。この第２周波数成分の振幅データおよび位相データを用いて作成した振幅画像および位相画像を図２６および図２７にそれぞれ示す。
【０１００】
図２７において、内膜部分が、子宮の体部から頚部に向かって徐々に暗くなっていることが観察できる。したがって、運動方向は、排卵期とは逆で、体部から頚部に向かう方向であると判断できる。月経期の子宮は、月経血の排出のため、体部から頚部に向かって運動していると考えられており、この実験では、そのことが確かめられた。
内膜上部の運動の速度を図２７中の点９と点８の画素の位相から求め、内膜下部の運動の速度を図２７中の点ｈと点ｅの位相から求めた。求めた速度を表２に示す。
【０１０１】
【表２】
【０１０２】
位相勾配を求める方向を、図２７の点８から点９向かう方向に指定して、各画素ごとに隣接画素との輝度値の差（位相勾配）を求め、正の位相勾配（点８から点９に向かって位相が増加する場合に対応）に対して青色の色データを対応付け、負の位相勾配（点９から点８に向かって位相が増加する場合に対応）に対して赤色の色データを対応付け、これらの色データを画素値とする位相勾配分布画像を作成すると、図２８(a)に示すカラー画像が得られる（(b)は(a)のカラー画像中の赤色部分を強調表示した画像。右下部の矢印は、左下に向かうものが青色矢印であり、右上に向かうものが赤色矢印である。）。すなわち、子宮内膜に対応する部分が赤色に着色された画像が得られ、子宮体部から頚部に向かう蠕動運動が明瞭に可視化されている。
【０１０３】
このカラー画像を原データに重ね合わせることにより、図２９に示すオーバレイ画像が得られる（ (a)はカラー画像、(b)は(a)のカラー画像の赤色領域を白抜きで強調表示した画像）。
［実験３（黄体期の画像］
図３０に示す対象画像（６４枚の子宮ＭＲ画像。図３０には代表的な１枚のみを示す。）は、或る被験者の黄体期の子宮の画像である。この画像の内膜の運動を解析し、黄体期の子宮内膜が、排卵期、月経期とはどのように異なるかを考察した。
【０１０４】
図３１(a)に、当該対象画像の最大振幅周波数分布画像を示す。参考のため、図３１(b)には図３１(a)のカラー画像中の第１周波数の色部分を強調した画像を示し、図３１(c)には図３１(a)のカラー画像中の第２周波数の色部分を強調した画像を示す。
図３１(a)の画像から、主要な運動の周波数成分（各部で最大振幅をとる周波数成分）は、第１周波数成分であると判断し、その振幅データおよび位相データに基づいて、振幅画像（図３２）および位相画像（図３３）をそれぞれ作成した。
【０１０５】
図３２の振幅画像には、内膜部分に高輝度な場所が明確には存在していない。また、図３３から、内膜部分の輝度値が一定であるため、運動による位相の変化は存在していないことが分かる。
図３４は、図３０中の４点ｑ，ｒ，ｓ，ｔの原データ（輝度値の時間変化）であり、図３５は、フーリエ変換によって得られた点ｑ，ｒ，ｓ，ｔの周波数スペクトル（各周波数成分の振幅データ）である。図３４から、時間軸方向の周期的な輝度値の変化がないことがわかる。また、図３５の振幅データから、最大振幅をとる周波数成分が明確ではなく、どの周波数成分の振幅の値も小さいことが分かる。
【０１０６】
これらの結果より、黄体期では、排卵期や月経期のような内膜運動は行われていないと判断できる。黄体期(分泌期) には、妊卵の着床や初期の妊卵の保持のため、内膜は蠕動運動を行っていないと考えられており、本実験による解析結果と一致する。
位相勾配を求める方向を、対称画像の右上角から左下角に向かう方向に指定して、各画素ごとに隣接画素との輝度値の差（位相勾配）を求め、正の位相勾配（右上角から左下角に向かって位相が増加する場合に対応）に対して青色の色データを対応付け、負の位相勾配（左下角から右上角に向かって位相が増加する場合に対応）に対して赤色の色データを対応付け、これらの色データを画素値とする位相勾配分布画像を作成すると、図３６に示す画像が得られる。この画像には、蠕動運動が現れていない。
［実験４（筋層の運動の解析）］
この実験では、子宮の内膜の外側を覆うように存在する子宮の筋層の運動に注目した。筋層は厚い筋肉から成り立っており、内膜と同様に周期的な運動をしていると考えられている。
【０１０７】
対象画像は図１４の場合（実験１）と同じである。したがって、最大振幅周波数分布画像は図１５(a)のとおりである。
この図１５(a)の最大振幅周波数分布画像からは、運動の周波が一定ではないと認められるため、筋層の運動を表す周波数成分を一意に決定することができない。そこで、いくつか周波数成分について、振幅画像と位相画像を作成し、図３７および図３８に振幅画像および位相画像をそれぞれ示す第６周波数成分が筋層に対応しているものと判断した（前述のとおり、第７周波数成分は内膜の蠕動に対応している）。
【０１０８】
図３８の位相画像より、運動方向は頚部から体部に向かう方向であると判断できる。筋層の運動は、内膜の運動よりも、画像中から手作業で捉えるのは一層困難であるが、本手法により解析が可能であることが確認できた。
排卵期の内膜の運動は頚部から体部に向かう方向の運動であったが、この実験により、筋層の運動方向もから体部に向かう方向であるという結果が得られた。つまり、内膜と筋層の運動方向は同じであるという結果が得られた。
【０１０９】
運動の速度を、図３８中の点ｕと点ｘから求めた。結果を表３に示す。
【０１１０】
【表３】
【０１１１】
次に、画像中で捉えられる筋層の運動から、運動の速度を手作業で求めた。画像列の１８枚目（図３９）から２４枚目（図４０）の間で、黒い影が、筋層を中の点Ｓから点Ｔまで移動するのが認められた。この運動から目視で求めた速度を表４に示す。
【０１１２】
【表４】
【０１１３】
画像中から目視で捉えられる運動の速度と、本手法により捉えた運動の速度がほぼ同じ値となっているため、筋層の蠕動の解析にも本手法が有効であると言える．
［他の実施形態］
前述の実施形態では、子宮の蠕動に対応した周波数成分の選択のために、最大振幅周波数分布画像がディスプレイ１３に表示され、操作者の判断によって、振幅画像および位相画像を作成すべき周波数成分が選択されているが、このような周波数成分の選択を自動化してもよい。具体的には、たとえば、最大振幅周波数分布画像中における個々の周波数成分の画素の数（すなわち、面積）を計算し、この画素数（面積）が上位の（たとえば最大の）周波数成分を自動で選択する周波数成分選択処理部４６（図６参照）を設けるようにしてもよい（面積順位基準選択手段）。また、周波数成分選択処理部４６は、所定のしきい値以上の最大振幅データに対応した１つまたは複数個の周波数成分を、子宮の蠕動に対応する蓋然性が高い周波数成分として自動選択する処理を行うものであってもよい。
【０１１４】
また、前述の実施形態では、算出すべき位相勾配の方向を入力装置１１から操作者が指定する構成としているが、位相勾配の方向は、予め定めた一定の方向（右上角から左下角に向かう方向、右から左に向かう方向など）であってもよい。
その他、特許請求の範囲に記載された事項の範囲で種々の設計変更を施すことが可能である。
【図面の簡単な説明】
【０１１５】
【図１】子宮ＭＲ画像（矢状断面）の一例を示す。
【図２】子宮ＭＲ画像の例を示す。
【図３】時間軸方向の輝度値の変化をプロットした例を示す。
【図４（ａ）】最大振幅周波数分布画像（カラー画像）の一例を示す。
【図４（ｂ）】図４(a)のカラー画像の第１周波数の色部分を強調表示した画像を示す。
【図４（ｃ）】図４(a)のカラー画像の第２周波数の色部分を強調表示した画像を示す。
【図４（ｄ）】図４(a)のカラー画像の第６周波数の色部分を強調表示した画像を示す。
【図５】この発明の一実施形態に係る画像診断支援システムの基本的なハードウェア構成を説明するためのブロック図である。
【図６】前記画像処理装置として機能するコンピュータ本体の機能的な構成を説明するためのブロック図である。
【図７】子宮ＭＲ画像データを処理して子宮の蠕動を可視化するための処理の流れをまとめて示すフローチャートである。
【図８】シミュレーション実験のために作成した６４枚の時系列子宮ＭＲ画像を説明するための画像例である。
【図９】図８の位置ＡおよびＢの画素の輝度値の時間変化を示す。
【図１０】図８の子宮ＭＲ画像についての最大振幅周波数分布画像を示す。
【図１１】振幅画像を示す。
【図１２】位相画像を示す。
【図１３】位相勾配分布画像を示す。
【図１４】排卵期の子宮ＭＲ画像の例を示す。
【図１５（ａ）】図１４の画像の最大振幅周波数分布画像を示す。
【図１５（ｂ）】図１５(a)のカラー画像の第３周波数の色部分を強調表示した画像を示す。
【図１５（ｃ）】図１５(a)のカラー画像の第６周波数の色部分を強調表示した画像を示す。
【図１５（ｄ）】図１５(a)のカラー画像の第７周波数の色部分を強調表示した画像を示す。
【図１６】子宮内膜の蠕動に対応した振幅画像を示す。
【図１７】子宮内膜の蠕動に対応した位相画像を示す
【図１８】図１７の４つの点ａ，ｂ，ｃ，ｄの輝度値の時間軸方向の変化を示す（原データ）。
【図１９】図１７の４つの点ａ，ｂ，ｃ，ｄの輝度値の時間軸方向の変化を示す（直流成分と第７周波数成分に逆フーリエ変換を適用したもの）。
【図２０】図１４の子宮ＭＲ画像列の２８枚目の子宮ＭＲ画像を示す。
【図２１】図１４の子宮ＭＲ画像列の３７枚目の子宮ＭＲ画像を示す。
【図２２】位相勾配分布画像を示す。
【図２３】オーバレイ画像を示す。
【図２４】月経期の子宮ＭＲ画像の例を示す。
【図２５】最大振幅周波数分布画像を示す。
【図２６】子宮内膜の蠕動に対応した振幅画像を示す。
【図２７】子宮内膜の蠕動に対応した位相画像を示す。
【図２８】位相勾配分布画像を示す。
【図２９】オーバレイ画像を示す。
【図３０】黄体期の子宮ＭＲ画像の例を示す。
【図３１】最大振幅周波数分布画像を示す。
【図３２】振幅画像を示す。
【図３３】位相画像を示す。
【図３４】図３０中の４点ｑ，ｒ，ｓ，ｔの原データ（輝度値の時間変化）を示す。
【図３５】図３０中の４点ｑ，ｒ，ｓ，ｔの周波数スペクトル（各周波数成分の振幅データ）である。
【図３６】位相勾配分布画像を示す。
【図３７】子宮筋層の蠕動に対応する振幅画像を示す。
【図３８】子宮筋層の蠕動に対応する位相画像を示す。
【図３９】図１４の子宮ＭＲ画像列の１８枚目の画像を示す。
【図４０】図１４の子宮ＭＲ画像列の２４枚目を示す。
【符号の説明】
【０１１６】
１０ コンピュータ本体
１１ 入力装置
１１Ａ キーボード
１１Ｂ マウス
１３ ディスプレイ
１４ プリンタ
１５ ＣＰＵ
１６ ＲＯＭ
１７ ＲＡＭ
１８ ハードディスクドライブ
１９ ＣＤ−ＲＯＭドライブ
２０ 入力制御部
２１ ディスプレイインタフェース
２２ プリンタインタフェース
２３ バス
３０ 画像処理演算部
３１ アフィン変換処理部
３２ トリミング処理部
３３ フーリエ変換処理部
３４ 最大振幅周波数成分抽出処理部
３５ 最大振幅周波数分布画像データ生成処理部
３６ 周波数成分選択入力受付け部
３７ 振幅画像データ生成処理部
３８ 位相画像データ生成処理部
３９ 位相勾配算出方向入力受付け部
４０ 位相勾配データ生成処理部
４１ 位相勾配分布画像データ生成処理部
４２ 画像オーバレイ処理部
４３ 速度測定点入力受付け部
４４ 位相差演算部
４５ 蠕動速度演算部
４６ 周波数成分選択処理部
５０ データ記憶部
５１ 原子宮ＭＲ画像データ記憶部
５２ アフィン変換後子宮ＭＲ画像データ記憶部
５３ 対象子宮ＭＲ画像データ記憶部
５４ フーリエ変換データ記憶部
５５ 最大振幅周波数成分データ記憶部
５６ 最大振幅周波数分布画像データ記憶部
５７ 振幅画像データ記憶部
５８ 位相画像データ記憶部
５９ 位相勾配データ記憶部
６０ 位相勾配分布画像データ記憶部
６１ オーバレイ画像データ記憶部
６２ 位相差データ記憶部
６３ 蠕動速度データ記憶部
７０ 画像表示制御部
【特許請求の範囲】
【請求項１】
子宮の蠕動運動を解析して可視化するための画像処理装置であって、
子宮の蠕動の周期よりも短い周期で時系列に従って撮影した複数枚の子宮ＭＲ（磁気共鳴）画像に対応した複数枚の子宮ＭＲ画像データを記憶するための画像データ記憶手段と、
この画像データ記憶手段に記憶された複数枚の子宮ＭＲ画像データに対して時間軸方向のフーリエ変換を行い、子宮ＭＲ画像中の各画素について複数の周波数成分の位相データを生成するフーリエ変換手段と、
このフーリエ変換手段によって生成された少なくとも１つの周波数成分の位相データを用いて、子宮ＭＲ画像各部における当該周波数成分の所定方向に関する位相勾配を算出し、その位相勾配を表す位相勾配データを生成する位相勾配データ生成手段と、
この位相勾配データ生成手段によって生成された位相勾配データに基づき、画像各部を位相勾配に応じて着色した位相勾配分布画像を表す位相勾配分布画像データを、子宮の蠕動を表す画像データとして生成する位相勾配分布画像データ生成手段とを含むことを特徴とする子宮蠕動可視化のための画像処理装置。
【請求項２】
前記画像データ記憶手段に記憶された複数枚の子宮ＭＲ画像データに対して、子宮画像部分の位置および形状の変動を補償するための幾何学変換を施し、この幾何学変換処理後の子宮ＭＲ画像データを生成する幾何学変換手段をさらに含み、
前記フーリエ変換手段は、前記幾何学変換手段による幾何学変換処理後の子宮ＭＲ画像データに対してフーリエ変換を実行するものであることを特徴とする請求項１記載の画像処理装置。
【請求項３】
前記フーリエ変換手段は、前記位相データのほかに、子宮ＭＲ画像の各画素についての複数の周波数成分の振幅データを生成するものであり、
このフーリエ変換手段が生成する振幅データに基づいて、少なくとも１つの周波数成分を選択する周波数成分選択手段をさらに含み、
前記位相勾配データ生成手段は、前記周波数成分選択手段によって選択された周波数成分に関して子宮ＭＲ画像各部の位相勾配データを生成するものであることを特徴とする請求項１または２記載の画像処理装置。
【請求項４】
前記周波数成分選択手段は、
子宮ＭＲ画像の各画素について、最大の振幅データに対応した最大振幅周波数成分を抽出する最大振幅周波数抽出手段と、
この最大振幅周波数抽出手段によって個々の画素について抽出された最大振幅周波数成分の子宮ＭＲ画像内での分布を表す最大振幅周波数分布画像データを生成する最大振幅周波数分布画像データ生成手段とを含むことを特徴とする請求項３記載の画像処理装置。
【請求項５】
前記周波数成分選択手段は、さらに、
操作者による周波数成分選択入力を受け付ける周波数成分選択入力受け付け手段を含むことを特徴とする請求項４記載の画像処理装置。
【請求項６】
前記周波数成分選択手段は、さらに、
前記最大振幅周波数分布画像データ生成手段が生成する最大振幅周波数分布画像データに基づき、当該最大振幅周波数分布画像データに対応した画像中に占める個々の周波数成分の面積の順位に基づいて、少なくとも１つの周波数成分を選択する面積順位基準選択手段を含むことを特徴とする請求項４記載の画像処理装置。
【請求項７】
前記周波数成分選択手段は、所定のしきい値以上の最大振幅データに対応した周波数成分を選択するものであることを特徴とする請求項３記載の画像処理装置。
【請求項８】
前記フーリエ変換手段は、前記位相データのほかに、子宮ＭＲ画像の各画素についての複数の周波数成分の振幅データを生成するものであり、
前記フーリエ変換手段によって生成された少なくとも１つの周波数成分の振幅データを用いて、子宮ＭＲ画像各部の当該周波数成分の振幅分布を表す振幅画像を表す振幅画像データを生成する振幅画像データ生成手段をさらに含むことを特徴とする請求項１ないし７のいずれかに記載の画像処理装置。
【請求項９】
前記フーリエ変換手段によって生成された少なくとも１つの周波数成分の位相データを用いて、子宮ＭＲ画像各部の当該周波数成分の位相分布を表す位相画像を表す位相画像データを生成する位相画像データ生成手段をさらに含むことを特徴とする請求項１ないし８のいずれかに記載の画像処理装置。
【請求項１０】
前記画像データ記憶手段に記憶された子宮ＭＲ画像データに前記位相勾配分布画像データ生成手段によって生成される位相勾配分布画像データをオーバレイして、オーバレイ画像データを生成する画像オーバレイ処理手段をさらに含むことを特徴とする請求項１ないし９のいずれかに記載の画像処理装置。
【請求項１１】
前記フーリエ変換手段によって生成された１つの周波数成分の位相データのうち、少なくとも２つの離間した画素の位相データに基づいて、当該２つの画素間の位相差を演算する位相差演算手段と、
この位相差演算手段によって演算された位相差に基づいて、子宮の蠕動速度を計算する蠕動速度演算手段をさらに含むことを特徴とする請求項１ないし１０のいずれかに記載の画像処理装置。
【請求項１２】
操作者による画像中の２点の指示入力を受け付ける手段をさらに含み、
前記位相差演算手段は、前記受け付けられた指示入力による２点の画素の位相データに基づいて、前記２点の画素間の位相差を演算するものであることを特徴とする請求項１１記載の画像処理装置。
【請求項１３】
請求項１ないし１２のいずれかに記載の画像処理装置と、
この画像処理装置に接続され、前記位相勾配分布画像を表示する表示装置とを含むことを特徴とする画像診断支援システム。
【請求項１４】
操作者によって操作可能であり、前記画像処理装置に接続されて、この画像処理装置に対して操作者の指示入力を入力するための入力装置をさらに含むことを特徴とする請求項１３記載の画像診断支援システム。
【請求項１５】
子宮の蠕動運動を解析して可視化するための画像処理方法であって、
子宮の蠕動の周期よりも短い周期で時系列に従って撮影した複数枚の子宮ＭＲ（磁気共鳴）画像に対応した複数枚の子宮ＭＲ画像データを準備するステップと、
この複数枚の子宮ＭＲ画像データに対して時間軸方向のフーリエ変換を行い、子宮ＭＲ画像中の各画素について複数の周波数成分の位相データを生成するフーリエ変換ステップと、
このフーリエ変換ステップによって生成された少なくとも１つの周波数成分の位相データを用いて、子宮ＭＲ画像各部における当該周波数成分の所定方向に関する位相勾配を算出し、その位相勾配を表す位相勾配データを生成する位相勾配データ生成ステップと、
この位相勾配データ生成ステップによって生成された位相勾配データに基づき、画像各部を位相勾配に応じて着色した位相勾配分布画像を表す位相勾配分布画像データを、子宮の蠕動を表す画像データとして生成する位相勾配分布画像データ生成ステップとを含むことを特徴とする子宮蠕動可視化のための画像処理方法。
【請求項１６】
前記フーリエ変換ステップによって生成された１つの周波数成分の位相データのうち、少なくとも２つの離間した画素の位相データに基づいて、当該２つの画素間の位相差を演算する位相差演算ステップと、
この位相差演算ステップによって演算された位相差に基づいて、子宮の蠕動速度を計算する蠕動速度演算ステップとをさらに含むことを特徴とする請求項１５記載の画像処理方法。
【請求項１７】
子宮の蠕動の周期よりも短い周期で時系列に従って撮影した複数枚の子宮ＭＲ（磁気共鳴）画像に対応した複数枚の子宮ＭＲ画像データを処理し、子宮の蠕動運動を解析して可視化するための画像処理装置としてコンピュータを動作させるためのコンピュータプログラムであって、
前記コンピュータを、
前記複数枚の子宮ＭＲ画像データに対して時間軸方向のフーリエ変換を行い、子宮ＭＲ画像中の各画素について複数の周波数成分の位相データを生成するフーリエ変換手段、
このフーリエ変換手段によって生成された少なくとも１つの周波数成分の位相データを用いて、子宮ＭＲ画像各部における当該周波数成分の所定方向に関する位相勾配を算出し、その位相勾配を表す位相勾配データを生成する位相勾配データ生成手段、および
この位相勾配データ生成手段によって生成された位相勾配データに基づき、画像各部を位相勾配に応じて着色した位相勾配分布画像を表す位相勾配分布画像データを、子宮の蠕動を表す画像データとして生成する位相勾配分布画像データ生成手段
として機能させることを特徴とする子宮蠕動可視化のためのコンピュータプログラム。
【請求項１８】
前記コンピュータを、さらに、
前記フーリエ変換手段によって生成された１つの周波数成分の位相データのうち、少なくとも２つの離間した画素の位相データに基づいて、当該２つの画素間の位相差を演算する位相差演算手段、および
この位相差演算手段によって演算された位相差に基づいて、子宮の蠕動速度を計算する蠕動速度演算手段
として機能させることを特徴とする請求項１７記載のコンピュータプログラム。
【請求項１９】
請求項１７または１８記載のコンピュータプログラムを記録しており、コンピュータによる読み取りが可能な記録媒体。
【請求項１】
子宮の蠕動運動を解析して可視化するための画像処理装置であって、
子宮の蠕動の周期よりも短い周期で時系列に従って撮影した複数枚の子宮ＭＲ（磁気共鳴）画像に対応した複数枚の子宮ＭＲ画像データを記憶するための画像データ記憶手段と、
この画像データ記憶手段に記憶された複数枚の子宮ＭＲ画像データに対して時間軸方向のフーリエ変換を行い、子宮ＭＲ画像中の各画素について複数の周波数成分の位相データを生成するフーリエ変換手段と、
このフーリエ変換手段によって生成された少なくとも１つの周波数成分の位相データを用いて、子宮ＭＲ画像各部における当該周波数成分の所定方向に関する位相勾配を算出し、その位相勾配を表す位相勾配データを生成する位相勾配データ生成手段と、
この位相勾配データ生成手段によって生成された位相勾配データに基づき、画像各部を位相勾配に応じて着色した位相勾配分布画像を表す位相勾配分布画像データを、子宮の蠕動を表す画像データとして生成する位相勾配分布画像データ生成手段とを含むことを特徴とする子宮蠕動可視化のための画像処理装置。
【請求項２】
前記画像データ記憶手段に記憶された複数枚の子宮ＭＲ画像データに対して、子宮画像部分の位置および形状の変動を補償するための幾何学変換を施し、この幾何学変換処理後の子宮ＭＲ画像データを生成する幾何学変換手段をさらに含み、
前記フーリエ変換手段は、前記幾何学変換手段による幾何学変換処理後の子宮ＭＲ画像データに対してフーリエ変換を実行するものであることを特徴とする請求項１記載の画像処理装置。
【請求項３】
前記フーリエ変換手段は、前記位相データのほかに、子宮ＭＲ画像の各画素についての複数の周波数成分の振幅データを生成するものであり、
このフーリエ変換手段が生成する振幅データに基づいて、少なくとも１つの周波数成分を選択する周波数成分選択手段をさらに含み、
前記位相勾配データ生成手段は、前記周波数成分選択手段によって選択された周波数成分に関して子宮ＭＲ画像各部の位相勾配データを生成するものであることを特徴とする請求項１または２記載の画像処理装置。
【請求項４】
前記周波数成分選択手段は、
子宮ＭＲ画像の各画素について、最大の振幅データに対応した最大振幅周波数成分を抽出する最大振幅周波数抽出手段と、
この最大振幅周波数抽出手段によって個々の画素について抽出された最大振幅周波数成分の子宮ＭＲ画像内での分布を表す最大振幅周波数分布画像データを生成する最大振幅周波数分布画像データ生成手段とを含むことを特徴とする請求項３記載の画像処理装置。
【請求項５】
前記周波数成分選択手段は、さらに、
操作者による周波数成分選択入力を受け付ける周波数成分選択入力受け付け手段を含むことを特徴とする請求項４記載の画像処理装置。
【請求項６】
前記周波数成分選択手段は、さらに、
前記最大振幅周波数分布画像データ生成手段が生成する最大振幅周波数分布画像データに基づき、当該最大振幅周波数分布画像データに対応した画像中に占める個々の周波数成分の面積の順位に基づいて、少なくとも１つの周波数成分を選択する面積順位基準選択手段を含むことを特徴とする請求項４記載の画像処理装置。
【請求項７】
前記周波数成分選択手段は、所定のしきい値以上の最大振幅データに対応した周波数成分を選択するものであることを特徴とする請求項３記載の画像処理装置。
【請求項８】
前記フーリエ変換手段は、前記位相データのほかに、子宮ＭＲ画像の各画素についての複数の周波数成分の振幅データを生成するものであり、
前記フーリエ変換手段によって生成された少なくとも１つの周波数成分の振幅データを用いて、子宮ＭＲ画像各部の当該周波数成分の振幅分布を表す振幅画像を表す振幅画像データを生成する振幅画像データ生成手段をさらに含むことを特徴とする請求項１ないし７のいずれかに記載の画像処理装置。
【請求項９】
前記フーリエ変換手段によって生成された少なくとも１つの周波数成分の位相データを用いて、子宮ＭＲ画像各部の当該周波数成分の位相分布を表す位相画像を表す位相画像データを生成する位相画像データ生成手段をさらに含むことを特徴とする請求項１ないし８のいずれかに記載の画像処理装置。
【請求項１０】
前記画像データ記憶手段に記憶された子宮ＭＲ画像データに前記位相勾配分布画像データ生成手段によって生成される位相勾配分布画像データをオーバレイして、オーバレイ画像データを生成する画像オーバレイ処理手段をさらに含むことを特徴とする請求項１ないし９のいずれかに記載の画像処理装置。
【請求項１１】
前記フーリエ変換手段によって生成された１つの周波数成分の位相データのうち、少なくとも２つの離間した画素の位相データに基づいて、当該２つの画素間の位相差を演算する位相差演算手段と、
この位相差演算手段によって演算された位相差に基づいて、子宮の蠕動速度を計算する蠕動速度演算手段をさらに含むことを特徴とする請求項１ないし１０のいずれかに記載の画像処理装置。
【請求項１２】
操作者による画像中の２点の指示入力を受け付ける手段をさらに含み、
前記位相差演算手段は、前記受け付けられた指示入力による２点の画素の位相データに基づいて、前記２点の画素間の位相差を演算するものであることを特徴とする請求項１１記載の画像処理装置。
【請求項１３】
請求項１ないし１２のいずれかに記載の画像処理装置と、
この画像処理装置に接続され、前記位相勾配分布画像を表示する表示装置とを含むことを特徴とする画像診断支援システム。
【請求項１４】
操作者によって操作可能であり、前記画像処理装置に接続されて、この画像処理装置に対して操作者の指示入力を入力するための入力装置をさらに含むことを特徴とする請求項１３記載の画像診断支援システム。
【請求項１５】
子宮の蠕動運動を解析して可視化するための画像処理方法であって、
子宮の蠕動の周期よりも短い周期で時系列に従って撮影した複数枚の子宮ＭＲ（磁気共鳴）画像に対応した複数枚の子宮ＭＲ画像データを準備するステップと、
この複数枚の子宮ＭＲ画像データに対して時間軸方向のフーリエ変換を行い、子宮ＭＲ画像中の各画素について複数の周波数成分の位相データを生成するフーリエ変換ステップと、
このフーリエ変換ステップによって生成された少なくとも１つの周波数成分の位相データを用いて、子宮ＭＲ画像各部における当該周波数成分の所定方向に関する位相勾配を算出し、その位相勾配を表す位相勾配データを生成する位相勾配データ生成ステップと、
この位相勾配データ生成ステップによって生成された位相勾配データに基づき、画像各部を位相勾配に応じて着色した位相勾配分布画像を表す位相勾配分布画像データを、子宮の蠕動を表す画像データとして生成する位相勾配分布画像データ生成ステップとを含むことを特徴とする子宮蠕動可視化のための画像処理方法。
【請求項１６】
前記フーリエ変換ステップによって生成された１つの周波数成分の位相データのうち、少なくとも２つの離間した画素の位相データに基づいて、当該２つの画素間の位相差を演算する位相差演算ステップと、
この位相差演算ステップによって演算された位相差に基づいて、子宮の蠕動速度を計算する蠕動速度演算ステップとをさらに含むことを特徴とする請求項１５記載の画像処理方法。
【請求項１７】
子宮の蠕動の周期よりも短い周期で時系列に従って撮影した複数枚の子宮ＭＲ（磁気共鳴）画像に対応した複数枚の子宮ＭＲ画像データを処理し、子宮の蠕動運動を解析して可視化するための画像処理装置としてコンピュータを動作させるためのコンピュータプログラムであって、
前記コンピュータを、
前記複数枚の子宮ＭＲ画像データに対して時間軸方向のフーリエ変換を行い、子宮ＭＲ画像中の各画素について複数の周波数成分の位相データを生成するフーリエ変換手段、
このフーリエ変換手段によって生成された少なくとも１つの周波数成分の位相データを用いて、子宮ＭＲ画像各部における当該周波数成分の所定方向に関する位相勾配を算出し、その位相勾配を表す位相勾配データを生成する位相勾配データ生成手段、および
この位相勾配データ生成手段によって生成された位相勾配データに基づき、画像各部を位相勾配に応じて着色した位相勾配分布画像を表す位相勾配分布画像データを、子宮の蠕動を表す画像データとして生成する位相勾配分布画像データ生成手段
として機能させることを特徴とする子宮蠕動可視化のためのコンピュータプログラム。
【請求項１８】
前記コンピュータを、さらに、
前記フーリエ変換手段によって生成された１つの周波数成分の位相データのうち、少なくとも２つの離間した画素の位相データに基づいて、当該２つの画素間の位相差を演算する位相差演算手段、および
この位相差演算手段によって演算された位相差に基づいて、子宮の蠕動速度を計算する蠕動速度演算手段
として機能させることを特徴とする請求項１７記載のコンピュータプログラム。
【請求項１９】
請求項１７または１８記載のコンピュータプログラムを記録しており、コンピュータによる読み取りが可能な記録媒体。
【図５】
【図６】
【図７】
【図９】
【図１８】
【図１９】
【図３４】
【図３５】
【図１】
【図２】
【図３】
【図４（ａ）】
【図４（ｂ）】
【図４（ｃ）】
【図４（ｄ）】
【図８】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５（ａ）】
【図１５（ｂ）】
【図１５（ｃ）】
【図１５（ｄ）】
【図１６】
【図１７】
【図２０】
【図２１】
【図２２】
【図２３】
【図２４】
【図２５】
【図２６】
【図２７】
【図２８】
【図２９】
【図３０】
【図３１】
【図３２】
【図３３】
【図３６】
【図３７】
【図３８】
【図３９】
【図４０】
【図６】
【図７】
【図９】
【図１８】
【図１９】
【図３４】
【図３５】
【図１】
【図２】
【図３】
【図４（ａ）】
【図４（ｂ）】
【図４（ｃ）】
【図４（ｄ）】
【図８】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５（ａ）】
【図１５（ｂ）】
【図１５（ｃ）】
【図１５（ｄ）】
【図１６】
【図１７】
【図２０】
【図２１】
【図２２】
【図２３】
【図２４】
【図２５】
【図２６】
【図２７】
【図２８】
【図２９】
【図３０】
【図３１】
【図３２】
【図３３】
【図３６】
【図３７】
【図３８】
【図３９】
【図４０】
【公開番号】特開２００６−２６３２８２（Ｐ２００６−２６３２８２Ａ）
【公開日】平成１８年１０月５日（２００６．１０．５）
【国際特許分類】
【出願番号】特願２００５−８７９４９（Ｐ２００５−８７９４９）
【出願日】平成１７年３月２５日（２００５．３．２５）
【出願人】（８０３０００１０４）財団法人ひろしま産業振興機構 (70)
【Ｆターム（参考）】
【公開日】平成１８年１０月５日（２００６．１０．５）
【国際特許分類】
【出願日】平成１７年３月２５日（２００５．３．２５）
【出願人】（８０３０００１０４）財団法人ひろしま産業振興機構 (70)
【Ｆターム（参考）】
[ Back to top ]