JP2014128687A - 動態画像診断支援システム - Google Patents

Abstract

【課題】呼吸状態下で撮影された胸部動態画像から、肺血流機能が低下している箇所を医師が容易に識別することができる動態画像診断支援システムを提供する。
【解決手段】診断用コンソール３によれば、制御部３１は、撮影用コンソール２から送信された胸部の動態画像の各フレーム画像から心臓壁位置を取得して心臓壁位置の時間変化を算出するとともに、動態画像の各フレーム画像の肺野領域を複数の小ブロックに分割し、分割された各小ブロック内の画素信号値の時間変化を算出する。そして、それぞれの時間変化を時間軸方向のハイパスフィルタでフィルタリングし、各小ブロックについて、心臓壁位置の時間変化に対する小ブロック内の画素信号値の時間変化の位相遅れ時間を算出し、各小ブロックを位相遅れ時間に応じた輝度値で表示部３４に表示する。
【選択図】図１

Description

本発明は、動態画像診断支援システムに関する。

従来のフィルム／スクリーンや輝尽性蛍光体プレートを用いた放射線の静止画撮影及び診断に対し、ＦＰＤ（flat panel detector）等の半導体イメージセンサを利用して検査対象部位の動態画像を撮影し、診断に応用する試みがなされるようになってきている。具体的には、半導体イメージセンサの画像データの読取・消去の応答性の早さを利用し、半導体イメージセンサの読取・消去のタイミングと合わせて放射源からパルス状の放射線を連続照射し、１秒間に複数回の撮影を行って、検査対象部位の動態を撮影する。撮影により取得された一連の複数枚の画像を順次表示することにより、医師は検査対象部位の一連の動きを認識することが可能となる。

また、動態画像を見やすく表示するための各種技術も提案されている。例えば、特許文献１には、時系列的に連続して取得された動態画像を、時間方向を奥行きに取った「時空間３次元画像」として扱い、その横断像、矢状断像を作成し、観察可能に表示する技術が記載されている。

特開２００４−３０５４８７号公報

ところで、肺の診断においては、肺の機能（換気機能や肺血流機能）が低下している箇所がないかを観察することは重要である。しかし、医師が動態画像を観察して目視で機能の異常個所を認識することは難しい。特に、肺の呼吸運動や心臓の拍動には個人差があり、個人差を考慮しながら換気機能や肺血流機能の異常箇所を視認することは困難である。また、肺血流機能の観察用の動態撮影時には患者は息を止める必要がある。

本発明の課題は、呼吸状態下で撮影された胸部動態画像から、肺血流機能が低下している箇所を医師が容易に識別することができるようにすることである。

上記課題を解決するため、請求項１に記載の発明の動態画像診断支援システムは、
呼吸状態下の胸部を動態撮影して、前記胸部の動態を示す複数のフレーム画像を取得する撮影手段と、
前記複数のフレーム画像における画素信号値の時間変化を時間方向のハイパスフィルタでフィルタリングして血流に関する情報を生成する画像解析手段と、
を備える。

請求項２に記載の発明は、請求項１に記載の発明において、
前記画像解析手段は、前記複数のフレーム画像における画素信号値の時間変化として、前記複数のフレーム画像の各フレーム画像を複数の領域に分割し、各領域内の画素信号値の時間変化を用いる。

請求項３に記載の発明は、請求項１又は２に記載の発明において、
前記画像解析手段は、前記複数のフレーム画像における画素信号値の時間変化を時間方向のハイパスフィルタでフィルタリングして得た前記血流に関する情報と、前記複数のフレーム画像から心臓の拍動に応じて変化する指標値の時間変化を時間方向のハイパスフィルタでフィルタリングして得た情報と、に基づいて、前記指標値の時間変化に対する前記画素信号値の時間変化の位相遅れ時間を算出し表示部に表示する。

請求項４に記載の発明は、請求項１〜３の何れか一項に記載の発明において、
前記ハイパスフィルタのカットオフ周波数が０．７Ｈｚである。

本発明によれば、呼吸状態下で撮影された胸部動態画像から、肺血流機能が低下している箇所を医師が容易に識別することが可能となる。

本発明の実施の形態における動態画像診断支援システムの全体構成を示す図である。 図１の撮影用コンソールの制御部により実行される撮影制御処理を示すフローチャートである。 図１の診断用コンソールの制御部により実行される画像解析処理を示すフローチャートである。 図３のステップＳ１１において実行される換気解析処理を示すフローチャートである。 １つの呼吸サイクルにおいて撮影された複数の時間位相Ｔ（Ｔ＝ｔ０〜ｔ６）のフレーム画像を示す図である。 肺尖と横隔膜の垂直方向の距離の算出時に用いられる各部を説明するための図である。 肺野領域を分割した小ブロックの例を示す図である。 横隔膜の上下位置の時間変化、ある小ブロックの平均信号値の時間変化、及び位相遅れ時間を示す図である。 図１の診断用コンソールの表示部に表示される横隔膜の上下位置の時間変化に対する遅れ度合いマップ等の一例を示す図である。 図３のステップＳ１２において実行される肺血流解析処理を示すフローチャートである。 心臓壁位置の算出時に用いられる各部を説明するための図である。 心臓壁位置の基準位置の算出を説明するための図である。 心臓壁位置の時間変化、ある小ブロックの平均信号値の時間変化、及び位相遅れ時間を示す図である。 図１の診断用コンソールの表示部に表示される心臓壁位置の時間変化に対する遅れ度合いマップ等の一例を示す図である。 大動脈弓と肺動脈を模式的に示す図である。

以下、図面を参照して本発明に係る実施の形態を詳細に説明する。ただし、発明の範囲は、図示例に限定されない。

〔動態画像診断支援システム１００の構成 〕
まず、構成を説明する。
図１に、本実施の形態における動態画像診断支援システム１００の全体構成を示す。
図１に示すように、動態画像診断支援システム１００は、撮影装置１と、撮影用コンソール２とが通信ケーブル等により接続され、撮影用コンソール２と、診断用コンソール３とがＬＡＮ（Local Area Network）等の通信ネットワークＮＴを介して接続されて構成されている。動態画像診断支援システム１００を構成する各装置は、ＤＩＣＯＭ（Digital Image and Communications in Medicine）規格に準じており、各装置間の通信は、ＤＩＣＯＭに則って行われる。

〔撮影装置１の構成〕
撮影装置１は、例えば、呼吸運動に伴う肺の膨張及び収縮の形態変化、心臓の拍動等の、周期性（サイクル）を持つ胸部の動態を撮影する装置である。動態撮影は、人体の胸部に対し、Ｘ線等の放射線を連続照射して複数の画像を取得（即ち、連続撮影）することにより行う。この連続撮影により得られた一連の画像を動態画像と呼ぶ。また、動態画像を構成する複数の画像のそれぞれをフレーム画像と呼ぶ。
撮影装置１は、図１に示すように、放射線源１１、放射線照射制御装置１２、放射線検出部１３、読取制御装置１４、サイクル検出センサ１５、サイクル検出装置１６等を備えて構成されている。

放射線源１１は、放射線照射制御装置１２の制御に従って、被写体Ｍに対し放射線（Ｘ線）を照射する。
放射線照射制御装置１２は、撮影用コンソール２に接続されており、撮影用コンソール２から入力された放射線照射条件に基づいて放射線源１１を制御して放射線撮影を行う。撮影用コンソール２から入力される放射線照射条件は、例えば、連続照射時のパルスレート、パルス幅、パルス間隔、撮影開始／終了タイミング、Ｘ線管電流の値、Ｘ線管電圧の値、フィルタ種等である。パルスレートは、１秒あたりの放射線照射回数であり、後述するフレームレートと一致している。パルス幅は、放射線照射１回当たりの放射線照射時間である。パルス間隔は、連続撮影において、１回の放射線照射開始から次の放射線照射開始までの時間であり、後述するフレーム間隔と一致している。

放射線検出部１３は、ＦＰＤ等の半導体イメージセンサにより構成される。ＦＰＤは、例えば、ガラス基板等を有しており、基板上の所定位置に、放射線源１１から照射されて少なくとも被写体Ｍを透過した放射線をその強度に応じて検出し、検出した放射線を電気信号に変換して蓄積する複数の画素がマトリックス状に配列されている。各画素は、例えばＴＦＴ（Thin Film Transistor）等のスイッチング部により構成されている。

読取制御装置１４は、撮影用コンソール２に接続されている。読取制御装置１４は、撮影用コンソール２から入力された画像読取条件に基づいて放射線検出部１３の各画素のスイッチング部を制御して、当該各画素に蓄積された電気信号の読み取りをスイッチングしていき、放射線検出部１３に蓄積された電気信号を読み取ることにより、画像データを取得する。この画像データがフレーム画像である。そして、読取制御装置１４は、取得したフレーム画像を撮影用コンソール２に出力する。画像読取条件は、例えば、フレームレート、フレーム間隔、画素サイズ、画像サイズ（マトリックスサイズ）等である。フレームレートは、１秒あたりに取得するフレーム画像数であり、パルスレートと一致している。フレーム間隔は、連続撮影において、１回のフレーム画像の取得動作開始から次のフレーム画像の取得動作開始までの時間であり、パルス間隔と一致している。

ここで、放射線照射制御装置１２と読取制御装置１４は互いに接続され、互いに同期信号をやりとりして放射線照射動作と画像の読み取りの動作を同調させるようになっている。

サイクル検出センサ１５は、被写体Ｍの呼吸運動の状態を検出して検出情報をサイクル検出装置１６に出力する。サイクル検出センサ１５としては、例えば、呼吸モニタベルト、ＣＣＤ（Charge Coupled Device）カメラ、光学カメラ、スパイロメータ等を適用することができる。

サイクル検出装置１６は、サイクル検出センサ１５により入力された検出情報に基づいて、呼吸サイクル数、及び現在呼吸運動の１サイクル中のどの状態であるか（例えば、吸気、吸気から呼気の変換点、呼気、呼気から吸気の変換点のどの状態か）を検出し、検出結果（サイクル情報）を撮影用コンソール２の制御部２１に出力する。サイクル検出装置１６は、例えば、サイクル検出センサ１５（呼吸モニタベルト、ＣＣＤカメラ、光学カメラ、スパイロメータ等）により肺の状態が吸気から呼気への変換点であることを示す検出情報が入力されたタイミングを１サイクルの基点とし、次にこの状態が検出されるタイミングまでの間を１サイクルとして認識する。

〔撮影用コンソール２の構成〕
撮影用コンソール２は、放射線照射条件や画像読取条件を撮影装置１に出力して撮影装置１による放射線撮影及び放射線画像の読み取り動作を制御するとともに、撮影装置１により取得された動態画像を撮影技師によるポジショニングの確認や診断に適した画像であるか否かの確認用に表示する。
撮影用コンソール２は、図１に示すように、制御部２１、記憶部２２、操作部２３、表示部２４、通信部２５を備えて構成され、各部はバス２６により接続されている。

制御部２１は、ＣＰＵ（Central Processing Unit）、ＲＡＭ（Random Access Memory）等により構成される。制御部２１のＣＰＵは、操作部２３の操作に応じて、記憶部２２に記憶されているシステムプログラムや各種処理プログラムを読み出してＲＡＭ内に展開し、展開されたプログラムに従って後述する撮影制御処理を始めとする各種処理を実行し、撮影用コンソール２各部の動作や、撮影装置１の放射線照射動作及び読み取り動作を集中制御する。

記憶部２２は、不揮発性の半導体メモリやハードディスク等により構成される。記憶部２２は、制御部２１で実行される各種プログラムやプログラムにより処理の実行に必要なパラメータ、或いは処理結果等のデータを記憶する。例えば、記憶部２２は、図２に示す撮影制御処理を実行するための撮影制御処理プログラムを記憶している。また、記憶部２２は、検査対象部位に対応付けて放射線照射条件及び画像読取条件を記憶している。各種プログラムは、読取可能なプログラムコードの形態で格納され、制御部２１は、当該プログラムコードに従った動作を逐次実行する。

操作部２３は、カーソルキー、数字入力キー、及び各種機能キー等を備えたキーボードと、マウス等のポインティングデバイスを備えて構成され、キーボードに対するキー操作やマウス操作により入力された指示信号を制御部２１に出力する。また、操作部２３は、表示部２４の表示画面にタッチパネルを備えても良く、この場合、タッチパネルを介して入力された指示信号を制御部２１に出力する。

表示部２４は、ＬＣＤ（Liquid Crystal Display）やＣＲＴ（Cathode Ray Tube）等のモニタにより構成され、制御部２１から入力される表示信号の指示に従って、操作部２３からの入力指示やデータ等を表示する。

通信部２５は、ＬＡＮアダプタやモデムやＴＡ（Terminal Adapter）等を備え、通信ネットワークＮＴに接続された各装置との間のデータ送受信を制御する。

〔診断用コンソール３の構成〕
診断用コンソール３は、撮影用コンソール２から動態画像を取得し、取得した動態画像を表示して医師が読影診断するための動画像処理装置である。
診断用コンソール３は、図１に示すように、制御部３１、記憶部３２、操作部３３、表示部３４、通信部３５を備えて構成され、各部はバス３６により接続されている。

制御部３１は、ＣＰＵ、ＲＡＭ等により構成される。制御部３１のＣＰＵは、操作部３３の操作に応じて、記憶部３２に記憶されているシステムプログラムや、各種処理プログラムを読み出してＲＡＭ内に展開し、展開されたプログラムに従って、後述する画像解析処理を始めとする各種処理を実行し、診断用コンソール３各部の動作を集中制御する。制御部３１は、後述する画像解析処理を実行することにより画像解析手段を実現する。

記憶部３２は、不揮発性の半導体メモリやハードディスク等により構成される。記憶部３２は、制御部３１で画像解析処理を実行するための画像解析処理プログラムを始めとする各種プログラムやプログラムにより処理の実行に必要なパラメータ、或いは処理結果等のデータを記憶する。これらの各種プログラムは、読取可能なプログラムコードの形態で格納され、制御部３１は、当該プログラムコードに従った動作を逐次実行する。

操作部３３は、カーソルキー、数字入力キー、及び各種機能キー等を備えたキーボードと、マウス等のポインティングデバイスを備えて構成され、キーボードに対するキー操作やマウス操作により入力された指示信号を制御部３１に出力する。また、操作部３３は、表示部３４の表示画面にタッチパネルを備えても良く、この場合、タッチパネルを介して入力された指示信号を制御部３１に出力する。

表示手段としての表示部３４は、ＬＣＤやＣＲＴ等のモニタにより構成され、制御部３１から入力される表示信号の指示に従って、操作部３３からの入力指示やデータ等を表示する。

通信部３５は、ＬＡＮアダプタやモデムやＴＡ等を備え、通信ネットワークＮＴに接続された各装置との間のデータ送受信を制御する。

〔動態画像診断支援システム１００の動作〕
次に、上記動態画像診断支援システム１００における動作について説明する。

（撮影装置１、撮影用コンソール２の動作）
まず、撮影装置１、撮影用コンソール２による撮影動作について説明する。
図２に、撮影用コンソール２の制御部２１において実行される撮影制御処理を示す。撮影制御処理は、制御部２１と記憶部２２に記憶されている撮影制御処理プログラムとの協働により実行される。

まず、撮影技師により撮影用コンソール２の操作部２３が操作され、撮影対象（被写体Ｍ）の患者情報（患者の氏名、身長、体重、年齢、性別等）の入力が行われる（ステップＳ１）。

次いで、放射線照射条件が記憶部２２から読み出されて放射線照射制御装置１２に設定されるとともに、画像読取条件が記憶部２２から読み出されて読取制御装置１４に設定される（ステップＳ２）。フレームレート（パルスレート）としては、５フレーム／秒以上が好ましい。通常、肺野内の換気、肺血流の遅延時間差は１秒以下であるため、位相遅れ時間（詳細後述）を精度良く複数の段階（少なくとも５段階以上）に分解して表現するには、５フレーム／秒以上を必要とするからである。

次いで、操作部２３の操作による放射線照射の指示が待機され、操作部２３により放射線照射指示が入力されると（ステップＳ３；ＹＥＳ）、サイクル検出装置１６にサイクル検出開始の指示が出力され、サイクル検出センサ１５及びサイクル検出装置１６による被写体Ｍの呼吸運動のサイクル検出が開始される（ステップＳ４）。サイクル検出装置１６により所定の状態（例えば、吸気から呼気への変換点）であることが検出されると、放射線照射制御装置１２及び読取制御装置１４に撮影開始指示が出力され、動態撮影が開始される（ステップＳ５）。即ち、放射線照射制御装置１２に設定されたパルス間隔で放射線源１１により放射線が照射され、放射線検出部１３によりフレーム画像が取得される。サイクル検出装置１６により予め定められた動態サイクル数が検出されると、制御部２１により放射線照射制御装置１２及び読取制御装置１４に撮影終了の指示が出力され、撮影動作が停止される。

撮影により取得されたフレーム画像は順次撮影用コンソール２に入力され、撮影順を示す番号と対応付けて記憶部２２に記憶されるとともに（ステップＳ６）、表示部２４に表示される（ステップＳ７）。撮影技師は、表示された動態画像によりポジショニング等を確認し、撮影により診断に適した画像が取得された（撮影ＯＫ）か、再撮影が必要（撮影ＮＧ）か、を判断する。そして、操作部２３を操作して、判断結果を入力する。

操作部２３の所定の操作により撮影ＯＫを示す判断結果が入力されると（ステップＳ８；ＹＥＳ）、動態撮影で取得された一連のフレーム画像のそれぞれに、動態画像を識別するための識別ＩＤや、患者情報、検査対象部位、放射線照射条件、画像読取条件、撮影順を示す番号、サイクル情報等の情報が付帯され（例えば、ＤＩＣＯＭ形式で画像データのヘッダ領域に書き込まれ）、通信部２５を介して診断用コンソール３に送信される（ステップＳ９）。そして、本処理は終了する。一方、操作部２３の所定の操作により撮影ＮＧを示す判断結果が入力されると（ステップＳ８；ＮＯ）、記憶部２２に記憶された一連のフレーム画像が削除され（ステップＳ１０）、本処理は終了する。

（診断用コンソール３の動作）
次に、診断用コンソール３における動作について説明する。
診断用コンソール３においては、通信部３５を介して撮影用コンソール２から動態画像の一連のフレーム画像が受信されると、制御部３１と記憶部３２に記憶されている画像解析処理プログラムとの協働により図３に示す画像解析処理が実行される。

画像解析処理においては、まず、換気解析処理が実行される（ステップＳ１１）。

図４に、換気解析処理のフローを示す。換気解析処理は、制御部３１と記憶部３２に記憶されている換気解析処理プログラムとの協働により実行される。

まず、動態画像から横隔膜の上下位置の時間変化が算出される（ステップＳ１０１）。

横隔膜は、その上下運動によって肺の呼吸運動を促すものである。図５に、１つの呼吸サイクルにおいて撮影された複数の時間位相Ｔ（Ｔ＝ｔ０〜ｔ６）のフレーム画像を示す。図５に示すように、呼吸サイクルは、呼気期と吸気期により構成される。呼気期は、横隔膜が上がることによって肺から空気が排出され、図５に示すように肺野の領域が小さくなる。最大呼気位では、横隔膜の位置が最も高い状態となる。吸気期は、横隔膜が下がることにより肺に空気が取り込まれ、図５に示すように胸郭中の肺野の領域が大きくなる。最大吸気位では、横隔膜の位置が最も下がった状態となる。このように、胸部動態画像において、横隔膜の上下位置の変化は、肺の呼吸運動を示す指標値となる。

ここで、図５からもわかるように、肺尖の上下位置は呼吸運動の影響をほとんど受けず、その位置がほとんど変わらないため、肺尖と横隔膜の垂直方向の距離Ｄ（図６参照）は、横隔膜の上下方向の位置を表しているといえる。よって、動態画像の各フレーム画像の肺尖と横隔膜の垂直方向の距離Ｄを算出することで、横隔膜の上下位置の時間変化を取得することができる。

例えば、ステップＳ１０１においては、まず、各フレーム画像について以下の処理が行われる。
はじめに、図６に示す左肺野及び右肺野の肺野領域Ｒの抽出が行われる。肺野領域Ｒは、放射線（Ｘ線）の透過量が多いため、その周辺の領域に比べて画素信号値（濃度値。以下、信号値という）が高くなる。そこで、まず、フレーム画像の各画素の信号値から濃度ヒストグラムを作成し、判別分析法等によって閾値を求める。次いで、求められた閾値より高信号の領域を肺野領域Ｒの候補として抽出する。次いで、候補領域の境界付近でエッジ検出を行い、境界付近の小領域でエッジが最大となる点を境界に沿って抽出していく。そして、抽出されたエッジ点を多項式関数で近似して肺野領域Ｒの境界線を取得する。
肺野領域Ｒが抽出されると、左右それぞれの肺野領域Ｒから肺尖及び横隔膜の基準位置Ｐ１、Ｐ２が特定される。例えば、肺尖の基準位置Ｐ１を肺野領域Ｒの最上端の位置と予め定義しておき、肺野領域Ｒにおける垂直方向の最も上方にある位置を抽出することによって肺尖の基準位置Ｐ１を特定する。また、横隔膜の基準位置Ｐ２を横隔膜のカーブＣ（図６に点線で示す）の垂直方向の平均位置と予め定義しておき、肺野領域Ｒから横隔膜のカーブＣを抽出し、その垂直方向の平均位置を求め、求めた位置を横隔膜の基準位置Ｐ２として特定する。そして、特定された肺尖の基準位置Ｐ１と横隔膜の基準位置Ｐ２の垂直方向の位置（Ｙ座標）の間の距離Ｄを算出する。

動態画像の各フレーム画像から左肺野と右肺野それぞれの上記Ｐ１とＰ２の垂直方向の距離Ｄが算出されると、距離Ｄの時間変化が呼吸運動を示す指標値としての横隔膜の上下位置の時間変化として取得される。
図８に、動態画像から取得された横隔膜の上下位置の時間変化を点線で示す。図８において横軸は、動態撮影開始からの経過時間ｔを、縦軸は距離Ｄを示している。

次いで、横隔膜の上下位置の時間変化の周期Ｔ０及び角周波数ω０が算出される（ステップＳ１０２）。
ステップＳ１０２では、例えば、肺尖と横隔膜の距離Ｄの最大値（ピーク）から次の最大値（ピーク）までの時間が周期Ｔ０として算出される。また、角周波数ω０は、ω０＝２π／Ｔ０によって算出される。

次いで、各フレーム画像にワーピング処理（非線形歪変換処理）が施され、フレーム画像間で肺野領域Ｒの位置合わせが行われる（ステップＳ１０３）。ステップＳ１０３では、例えば、各フレーム画像において、同等の解剖学的特徴位置の複数の点（ランドマーク）を抽出し、１番目に撮影されたフレーム画像を基準画像として、基準画像からのランドマーク各点のシフト値を算出する。そして、算出されたランドマークの各点のシフト値に基づいて、全画素のシフト値を６次元多項式の近似処理を用いて算出し、算出されたシフト値分だけ各画素をシフトした画像を生成する。

次いで、各フレーム画像の肺野領域Ｒが複数の領域（小ブロックＡ１）に分割され、各小ブロックＡ１内の画素の平均信号値（濃度平均値）が算出される（ステップＳ１０４）。例えば、図７に示すように、肺野領域Ｒに外接する矩形領域が０．４〜４ｃｍ角の小ブロックＡ１に分割され、ブロック毎に画素の平均信号値が算出される。

次いで、各小ブロックＡ１ごとに、上記算出された平均信号値の時間変化が算出される（ステップＳ１０５）。
図８に、ある小ブロックＡ１の平均信号値の時間変化を実線で示す。図８において横軸は動態撮影開始からの経過時間ｔを、縦軸は、小ブロックＡ１内の画素の平均信号値を示している。放射線（Ｘ線）画像における肺野の信号値は肺野のＸ線の透過量によって定まり、Ｘ線の透過量は肺の膨張・収縮によって変わる。従って、肺野のある領域の信号値が変化しているということは肺の膨張・収縮に伴ってその領域の肺胞に含まれる空気の量が変化している、つまり換気量が変化していると考えることができ、各小ブロックＡ１の平均信号値の変化はその小ブロックＡ１における換気量の変化と相関関係がある。

次いで、各小ブロックＡ１ごとに、横隔膜の上下位置の時間変化（指標値の時間変化）に対する平均信号値の時間変化の位相遅れ時間が算出される（ステップＳ１０６）。
ステップＳ１０６における位相遅れ時間は、フーリエ級数展開による時間遅れ時間算出方法により算出することができる。
まず、指標値の時間変化の位相α１と、各小ブロックＡ１の平均信号値の時間変化の位相α２が算出される。位相α１は、下記の〔数１〕により算出される。
位相α２は、上記α１におけるｘ（ｔ）を時間ｔのときの小ブロックＡ１の平均信号値に置き換えることにより算出することができる。
位相α１、α２が算出されると、指標値の時間変化に対する各小ブロックＡ１の平均信号値の時間変化の位相遅れ時間αＴが算出される。αＴは、下記の〔数２〕により算出される。

または、ステップＳ１０６における位相遅れ時間は、相互相関により算出することができる。相互相関による位相遅れ時間の算出では、まず、各小ブロックＡ１の平均信号値の時間変化を時間方向のローパスフィルタ（例えば、カットオフ周波数０．５Ｈｚ）でフィルタリングする。次いで、指標値としての横隔膜の上下位置の時間変化と、フィルタリング後の各小ブロックＡ１の平均信号値の時間変化との相互相関値を算出する。そして、相互相関値が最大となる時間を各小ブロックＡ１の位相遅れ時間αＴとして算出する。

上記のステップＳ１０６の算出方法によれば、横隔膜の上下位置の時間変化と略同じ周波数の画素信号成分を抽出して位相遅れ時間を算出しているので、各小ブロックＡ１の平均信号値の時間変化から換気量についての情報のみを取り出し、肺血流の影響を取り除くことができ、精度良く換気量の位相遅れ時間を算出することが可能となる。

そして、各小ブロックＡ１の位相遅れ時間に応じた輝度値が算出され、基準画像の各小ブロックＡ１を上記算出された輝度値により表示した遅れ度合いマップＭ１が表示部３４ に表示される（ステップＳ１０７）。

また、記憶部３２から各小ブロックＡ１の位相遅れ時間の閾値が読み出され、ステップＳ１０６で算出された各小ブロックＡ１の位相遅れ時間が閾値を超えるか否かによって、各小ブロックＡ１の換気機能が異常であるか否かが判断される（ステップＳ１０８）。記憶部３２には、横隔膜の基準位置Ｐ２として特定された位置からの距離に応じた閾値が記憶されている。ステップＳ１０８では、各小ブロックＡ１について、基準位置Ｐ２からの距離が算出され、位相遅れ時間が算出された距離に応じた閾値を超えるか否かが判断される。
そして、上述の基準画像における異常と判断された小ブロックＡ１を他のブロックとは異なる色で表示した異常判定結果Ｍ２が遅れ度合いマップＭ１に並べて表示される（ステップＳ１０９）。

図９に、表示部３４に表示される遅れ度合いマップＭ１、異常判定結果Ｍ２の一例を示す。また、図９には、Ｍ１、Ｍ２と併せて標準遅れ度合いマップＭ０が表示されている。標準遅れ度合いマップは、正常な肺野における、横隔膜の上下位置の時間変化に対する各小ブロックＡ１の平均信号値の時間変化の位相遅れ時間を示したものである。この標準遅れ度合いマップＭ０は、正常な肺野を撮影した複数の動態画像に基づき作成されたものである。
図９の標準遅れ度合いマップＭ０及び遅れ度合いマップＭ１では、白（最高輝度値）は呼吸運動の指標値である横隔膜の上下位置の時間変化から遅れがないことを示し、所定ミリ秒（例えば、３３ミリ秒）遅れる毎に輝度値を１段階下げて表示される。

肺野では、横隔膜の上下運動によって呼吸運動が促され、換気が行われるため、図９の標準遅れ度合いマップＭ０に示されるように、横隔膜からの距離によって、各小ブロックＡ１の換気量の変化を示す平均信号値の時間変化の位相遅れ時間は異なる。しかし、その位相遅れ度合いは、換気機能が正常であれば、横隔膜からの距離に応じて略一定である。しかし、肺野領域内に換気機能が異常な箇所があると、その部分の位相遅れは大きくなる。そこで、図９に示すように、標準遅れ度合いマップＭ０と、被写体Ｍを動態撮影することにより得られた遅れ度合いマップＭ１とを並べて表示することで、被写体Ｍの肺野において局所的に換気機能が低下している箇所を医師が容易に識別することが可能となる。また、異常判定結果Ｍ２に、予め定められた閾値より位相遅れ時間の大きい小ブロックＡ１を異なる色で表示することで、局所的に換気機能が低下している異常箇所を医師がより容易に識別することが可能となる。遅れ度合いマップＭ１、異常判定結果Ｍ２を表示することで、例えば、慢性閉塞性肺疾患（ＣＯＰＤ）、間質性肺炎、気胸等がある場合に、医師は容易に認識することが可能となる。図９の遅れ度合いマップＭ１によれば、被写体Ｍは左肺野上部に気胸等による異常の疑いがあることが分かる。

操作部３３により、肺血流解析への移行が指示されると、処理は図３のステップＳ１２に移行し、肺血流解析処理が実行される（ステップＳ１２）。

図１０に、肺血流解析処理のフローを示す。当該処理は、制御部３１と記憶部３２に記憶されている肺血流解析処理プログラムとの協働により実現される。

まず、動態画像から左心室の心臓壁位置の時間変化が算出される（ステップＳ２０１）。
心臓壁の位置は、心臓の拍動によって変化するものであるため、胸部動態画像において、心臓壁の位置の時間変化は、心臓の拍動を示す指標値として用いることができる。

例えば、ステップＳ２０１においては、まず、各フレーム画像について以下の処理が行われる。
まず、図１１に示す左肺野及び右肺野の肺野領域Ｒの抽出、及び左肺野における肺尖の位置Ｐ１の特定が行われる。
次いで、心臓領域Ｈが認識される。心臓領域Ｈの認識においては、まず、肺野領域Ｒの外接矩形領域から探索領域を限定する。次いで、探索領域の各画素の信号値から濃度ヒストグラムを作成し、判別分析法等によって閾値を求め、閾値より低信号の領域を心臓の候補領域として抽出する。次いで、候補領域内でエッジ検出を行い、所定の大きさ以上の微分値の極大値を追跡することで、心臓領域の輪郭線を抽出する。このとき、背景若しくは心臓内部のエッジを追跡しないように、近似的心臓領域の形状をもとに輪郭エッジ点の探索領域を限定する。そして、抽出された心臓領域の輪郭線画像に対して、心臓輪郭線テンプレートでテンプレートマッチングを行い、相関値が最大となる位置でのテンプレート領域を心臓領域Ｈとして認識する。
次いで、左心室の心臓壁の基準位置Ｐ３が特定される。例えば、肺尖の基準位置Ｐ１から垂直方向が所定の距離にある心臓領域Ｈの最も左側のエッジ点を心臓壁の基準位置Ｐ３として特定する。なお、心臓領域Ｈの左側エッジの緩やかな湾曲部分が左心室に相当するので、上から下に向かって、垂直方向の座標（Ｙ座標）に対する心臓エッジの勾配（傾きの絶対値（図１２のΔＸ／ΔＹの絶対値）を求めていき、その勾配が、所定の閾値以下となってから初めて肺尖からの距離が所定の閾値以上となる点の位置を心臓壁の基準位置Ｐ３として特定してもよい。

動態画像の各フレーム画像から上記Ｐ３の位置が特定されると、各フレーム画像のＰ３の水平方向の位置（Ｐ３のＸ座標）に基づいて、心臓壁位置の時間変化が算出される。
図１３に、動態画像から算出された心臓壁位置の時間変化を点線で示す。図１３において横軸は、動態撮影開始からの経過時間ｔを、縦軸は心臓壁の基準位置Ｐ３の水平方向の位置を示している。

次いで、心臓壁位置の時間変化の周期Ｔ０及び角周波数ω０が算出される（ステップＳ２０２）。
ステップＳ２０２では、例えば、心臓壁位置の最大値（ピーク）から次の最大値（ピーク）までの時間が周期Ｔ０として算出される。また、角周波数ω０は、ω０＝２π／Ｔ０によって算出される。

次いで、各フレーム画像にワーピング処理（非線形歪変換処理）が施され、フレーム画像間で肺野領域Ｒの位置合わせが行われる（ステップＳ２０３）。ステップＳ２０３の処理は、ステップＳ１０３の処理と同じである。

次いで、各フレーム画像の肺野領域Ｒが小ブロックＡ１に分割され、各小ブロックＡ１内の画素の平均信号値（濃度平均値）が算出される（ステップＳ２０４）。例えば、肺野領域Ｒに外接する矩形領域が０．４〜４ｃｍ角の小ブロックＡ１に分割され、ブロック毎に画素の平均信号値が算出される。

次いで、各小ブロックＡ１ごとに、上記算出された平均信号値の時間変化が算出される（ステップＳ２０５）。
図１３に、ある小ブロックＡ１の平均信号値の時間変化を実線で示す。図１３において横軸は動態撮影開始からの経過時間ｔを、縦軸は、小ブロックＡ１内の画素の平均信号値を示している。放射線（Ｘ線）画像における肺野の信号値は、心臓の拍動によって生じる肺血流量によって変わる。従って、各小ブロックＡ１の平均信号値の変化はその小ブロックＡ１における肺血流量の変化と相関関係がある。

次いで、各小ブロックＡ１ごとに、心臓壁位置の時間変化（指標値の時間変化）に対する平均信号値の時間変化の位相遅れ時間が算出される（ステップＳ２０６）。
ステップＳ２０６における位相遅れ時間は、例えば、上述の〔数１〕〔数２〕に示すフーリエ級数展開による時間遅れ時間算出方法により算出することができる。または、相互相関により位相遅れ時間を算出することができる。相互相関による位相遅れ時間の算出では、まず、心臓壁位置の時間変化を時間方向のハイパスフィルタ（例えば、カットオフ周波数０．７Ｈｚ）でフィルタリングする。また、各小ブロックＡ１の平均信号値の時間変化を時間方向のハイパスフィルタ（例えば、カットオフ周波数０．７Ｈｚ）でフィルタリングする。次いで、フィルタリング後の心臓壁位置の時間変化と、フィルタリング後の各小ブロックＡ１の平均信号値の時間変化との相互相関値を算出する。そして、相互相関値が最大となる時間を各小ブロックＡ１の位相遅れ時間として算出する。

上記のステップＳ２０６の算出方法によれば、心臓壁位置の時間変化と略同じ周波数の画素信号成分を抽出して位相遅れ時間を算出しているので、各小ブロックＡ１の平均信号値の時間変化から肺血流量についての情報のみを取り出して、換気（呼吸）の影響を取り除くことができ、精度良く肺血流量の位相遅れ時間を算出することが可能となる。また、動態撮影時に肺血流解析用に息を止めることが不要となる。なお、息を止めた状態で動態撮影を行った場合は、呼吸の影響がないのでフィルタリングの必要はない。

そして、各小ブロックＡ１の位相遅れ時間に応じた輝度値が算出され、基準画像の各小ブロックＡ１を上記算出された輝度値により表示した遅れ度合いマップＭ１１が表示部３４に表示される（ステップＳ２０７）。

また、記憶部３２から各小ブロックＡ１の位相遅れ時間の閾値が読み出され、ステップＳ２０６で算出された各小ブロックＡ１の位相遅れ時間が閾値を超えるか否かによって、各小ブロックＡ１の肺血流機能が異常であるか否かが判断される（ステップＳ２０８）。記憶部３２には、心臓の中心部からの距離に応じた位相遅れ時間の閾値が記憶されている。ステップＳ２０８では、各小ブロックＡ１について、心臓の中心部（例えば、心臓領域Ｈの外接円の中心等）からの距離が算出され、位相遅れ時間が算出された距離に応じた閾値を超えるか否かが判断される。
そして、上述の基準画像における異常と判断された小ブロックＡ１を他のブロックとは異なる色で表示した異常判定結果Ｍ１２が遅れ度合いマップＭ１１に並べて表示される（ステップＳ２０９）。

図１４に、表示部３４に表示される遅れ度合いマップＭ１１、異常判定結果Ｍ１２の一例を示す。また、図１４には、Ｍ１１、Ｍ１２と併せて標準遅れ度合いマップＭ１０が表示されている。標準遅れ度合いマップＭ１０は、正常な肺野における、心臓の拍動の指標値である心臓壁位置の時間変化に対する各小ブロックＡ１の平均信号値の時間変化の位相遅れ時間を示したものである。この標準遅れ度合いマップＭ１０は、正常な肺野を撮影した複数の動態画像に基づき作成されたものである。
図１４の標準遅れ度合いマップＭ１０及び遅れ度合いマップＭ１１では、各小ブロックＡ１の位相遅れ時間を輝度で表示しているが、白（最高輝度値）は心臓の拍動の指標値である心臓壁の時間変化から遅れがないことを示し、所定ミリ秒（例えば、３３ミリ秒）遅れる毎に輝度値が１段階下げて表示される。

肺野では、心臓の拍動により心臓から血液が押し出され、この血液が肺野各部の血管を流れるため、図１４の標準遅れ度合いマップＭ１０に示されるように、心臓の中心部からの距離によって、各小ブロックＡ１の肺血流量を示す平均信号値の時間変化の位相遅れは異なる。しかし、その位相遅れ度合いは、肺血流機能が正常であれば、心臓からの距離に応じて略一定である。しかし、肺野領域内に肺血流機能が異常な箇所があると、その部分の位相遅れは大きくなる。そこで、図１４に示すように、標準遅れ度合いマップＭ１０と、患者を動態撮影することにより得られた遅れ度合いマップＭ１１とを並べて表示することで、局所的に肺血流機能が低下している箇所を医師が容易に識別することが可能となる。また、異常判定結果Ｍ１２に、予め定められた閾値より位相遅れ時間の大きい小ブロックＡ１を異なる色で表示することで、局所的に肺血流機能が低下している異常箇所を医師がより容易に識別することが可能となる。遅れ度合いマップＭ１１、異常判定結果Ｍ１２を表示することで、例えば、肺血栓塞栓症、肺梗塞症等の病変がある場合に、医師は容易に認識することが可能となる。図１４の遅れ度合いマップＭ２によれば、被写体Ｍは右肺野に肺血栓塞栓症等による異常の疑いがあることがわかる。

なお、上記の肺血流解析処理では、心臓の拍動を示す指標値として心臓壁位置を用い、心臓壁位置の時間変化に対する各小ブロックＡ１の平均信号値の時間変化の位相遅れ時間を輝度値等により表示することとしたが、心臓の拍動を示す指標値としては、これに限定されない。例えば、動態画像における大動脈弓又は肺動脈の部分の信号値を心臓の拍動を示す指標値として用い、大動脈弓又は肺動脈の部分の信号値に対する各小ブロックＡ１の時間変化の位相遅れ時間を輝度値等により表示することとしてもよい。

図１５に、胸部画像における大動脈弓、肺動脈を模式的に示す。
図１５に示すように、大動脈弓は、通常、胸部正面Ｘ線画像上では、左肺野縦隔側の中央よりやや上に見られ、形状が左向きに凸状に湾曲した陰影として現れる。胸部正面Ｘ線画像（即ち、各フレーム画像）からエッジを抽出した上で、大動脈弓陰影のテンプレートを用いてテンプレートマッチングすることで、その位置を認識することができる。この認識した大動脈弓陰影内にＲＯＩ（Region Of Interest）を設定し、ＲＯＩ内の平均信号値を心臓の拍動を示す指標値として用いる。
左肺動脈は、左肺野内の大動脈弓やや下に位置し、肺紋理と呼ばれる白い紋様の陰影に含まれ、その中でも太い動脈は縦隔部から若干左側に突き出て、下側に向かって走行する陰影として現れる。大動脈弓と同様、エッジを抽出した上で、この太い動脈の陰影のテンプレートを用いてテンプレートマッチングをすることで、その位置を認識することができる。認識した左肺動脈陰影内の極力縦隔に近い部分にＲＯＩを設定し、ＲＯＩ内の平均信号値を心臓の拍動を示す指標値として用いる。

大動脈弓は、心臓の左心室から全身に血液を送り出す動脈の本幹であり、肺動脈は、右心室から肺へ血液を送り出す動脈であり、何れの信号値も心臓の拍動を顕著に示すものである。肺血流量を観察する上では肺動脈上の信号値変化を抽出する方が好ましいが、肺動脈は大動脈弓に比べて陰影が淡いため識別し難く、また心臓に隠れる場合があるので、画像処理上は大動脈弓の信号値を心臓の拍動の指標値とすることが好ましい。

なお、画像解析処理の処理結果（各小ブロックＡ１の位相遅れ時間及び異常判定結果）は、動態画像と対応付けて記憶部３２に記憶しておくことが好ましい。次に動態画像を読影する際に、画像解析処理を行う必要がなくなるためである。

以上説明したように、本発明に係る診断用コンソール３によれば、制御部３１は、撮影用コンソール２から送信された胸部の動態画像の各フレーム画像から横隔膜の上下位置を取得して横隔膜の上下位置の時間変化を算出するとともに、動態画像の各フレーム画像の肺野領域を複数の小ブロックＡ１に分割し、分割された各小ブロックＡ１内の画素信号値の時間変化を算出し、各小ブロックＡ１について、横隔膜の上下位置の時間変化に対する小ブロック内の画素信号値の時間変化の位相遅れ時間を算出する。そして、各小ブロックＡ１を位相遅れ時間に応じた輝度値で表示部３４に表示する。

従って、肺野内において、肺の呼吸運動の指標となる横隔膜の上下位置の時間変化に対して、換気量を示す信号値の時間変化の位相遅れ時間がどのように分布しているかを表示するので、肺の換気機能が低下している箇所があった場合に医師が容易に識別することが可能となる。

また、制御部３１は、算出された各小ブロックＡ１の位相遅れ時間が予め定められた閾値を超えるか否かに基づいて各小ブロックＡ１の換気機能が異常であるか否かを判断し、判断結果を表示部３４に表示するので、肺の換気機能が低下している異常な箇所を医師がより容易に認識することが可能となる。また、被写体Ｍの呼吸運動を基準として、その基準に対し肺野内の各ブロックの信号値の変化にどの程度遅れがあるかによって各ブロックの換気機能の異常判断を行うので、呼吸運動の個人差を考慮して精度良く異常判断を行うことが可能となる。

位相遅れ時間の算出にあたっては、横隔膜の上下位置の時間変化と略同じ周波数の画素信号成分を抽出して位相遅れ時間を算出するので、動態画像から換気量についての情報のみを取り出し、肺血流の影響を取り除くことができ、精度良く位相遅れ時間を算出することが可能となる。

また、診断用コンソール３の制御部３１は、撮影用コンソール２から送信された胸部の動態画像の各フレーム画像から心臓壁位置を取得して心臓壁位置の時間変化を算出するとともに、動態画像の各フレーム画像の肺野領域を複数の小ブロックＡ１に分割し、分割された各小ブロックＡ１内の画素信号値の時間変化を算出し、各小ブロックＡ１について、心臓壁位置の時間変化に対する小ブロック内の画素信号値の時間変化の位相遅れ時間を算出する。そして、各小ブロックＡ１を位相遅れ時間に応じた輝度値で表示部３４に表示する。

従って、肺野内において、心臓の拍動の指標となる心臓壁位置の時間変化に対して、肺血流量を示す信号値の時間変化の位相遅れ時間がどのように分布しているかを表示するので、肺血流機能が低下している箇所があった場合に医師が容易に識別することが可能となる。

また、制御部３１は、算出された各小ブロックＡ１の位相遅れ時間が予め定められた閾値を超えるか否かに基づいて各小ブロックＡ１の肺血流機能が異常であるか否かを判断し、判断結果を表示部３４に表示するので、肺血流機能が低下している異常な箇所を医師がより容易に認識することが可能となる。また、被写体Ｍの心臓の拍動を基準として、その基準に対し肺野内の各ブロックの信号値の変化にどの程度遅れがあるかによって各ブロックの肺血流機能の異常判断を行うので、心臓の拍動の個人差を考慮して精度良く異常判断を行うことが可能となる。

位相遅れ時間の算出にあたっては、心臓壁位置の時間変化と略同じ周波数の画素信号成分を抽出して位相遅れ時間を算出するので、動態画像から肺血流量についての情報のみを取り出し、換気の影響を取り除くことができ、精度良く位相遅れ時間を算出することが可能となる。また、撮影時に被写体Ｍが息を止める必要がなくなる。

また、画像解析処理においては、呼吸時に撮影された胸部の動態画像から、横隔膜の上下位置の時間変化に対する各小ブロックＡ１の平均信号値の時間変化の位相遅れ時間、及び心臓壁位置の時間変化に対する各小ブロックＡ１の平均信号値の時間変化の位相遅れ時間をそれぞれ算出し、各小ブロックＡ１について算出されたそれぞれの位相遅れ時間をそれぞれ表示部３４に表示するので、一回の撮影で、換気機能と肺血流機能の双方の解析を行うことが可能となる。

なお、上述した本実施の形態における記述は、本発明に係る好適な動態画像診断支援システムの一例であり、これに限定されるものではない。
例えば、上記実施の形態においては、標準遅れ度合いマップＭ０、Ｍ１０、遅れ度合いマップＭ１、Ｍ１１を位相遅れ時間に応じた輝度値により表示することとしたが、位相遅れ時間に応じた色で表示することとしてもよい。

また、上記実施の形態においては、換気機能と肺血流機能の双方を解析する場合を例に取り説明したが、操作部３３の操作により、解析対象とする機能を選択できるようにしてもよい。

また、例えば、上記の説明では、本発明に係るプログラムのコンピュータ読み取り可能な媒体としてハードディスクや半導体の不揮発性メモリ等を使用した例を開示したが、この例に限定されない。その他のコンピュータ読み取り可能な媒体として、ＣＤ-ＲＯＭ等の可搬型記録媒体を適用することが可能である。また、本発明に係るプログラムのデータを通信回線を介して提供する媒体として、キャリアウエーブ(搬送波)も適用される。

その他、動態画像診断支援システム１００を構成する各装置の細部構成及び細部動作に関しても、本発明の趣旨を逸脱することのない範囲で適宜変更可能である。

１００ 動態画像診断支援システム
１ 撮影装置
１１ 放射線源
１２ 放射線照射制御装置
１３ 放射線検出部
１４ 読取制御装置
１５ サイクル検出センサ
１６ サイクル検出装置
２ 撮影用コンソール
２１ 制御部
２２ 記憶部
２３ 操作部
２４ 表示部
２５ 通信部
２６ バス
３ 診断用コンソール
３１ 制御部
３２ 記憶部
３３ 操作部
３４ 表示部
３５ 通信部
３６ バス

Claims (4)

1. 呼吸状態下の胸部を動態撮影して、前記胸部の動態を示す複数のフレーム画像を取得する撮影手段と、
   前記複数のフレーム画像における画素信号値の時間変化を時間方向のハイパスフィルタでフィルタリングして血流に関する情報を生成する画像解析手段と、
   を備える動態画像診断支援システム。
2. 前記画像解析手段は、前記複数のフレーム画像における画素信号値の時間変化として、前記複数のフレーム画像の各フレーム画像を複数の領域に分割し、各領域内の画素信号値の時間変化を用いる請求項１に記載の動態画像診断支援システム。
3. 前記画像解析手段は、前記複数のフレーム画像における画素信号値の時間変化を時間方向のハイパスフィルタでフィルタリングして得た前記血流に関する情報と、前記複数のフレーム画像から心臓の拍動に応じて変化する指標値の時間変化を時間方向のハイパスフィルタでフィルタリングして得た情報と、に基づいて、前記指標値の時間変化に対する前記画素信号値の時間変化の位相遅れ時間を算出し表示部に表示する請求項１又は２に記載の動態画像診断支援システム。
4. 前記ハイパスフィルタのカットオフ周波数が０．７Ｈｚである請求項１〜３の何れか一項に記載の動態画像診断支援システム。