(57)【特許請求の範囲】
【請求項１】
検査対象部位に放射線を照射して動態撮影することにより得られた動態画像に対して２軸以上の特徴量空間でそれぞれの軸の特徴量でフィルタリングを施すためのパラメーターを設定するパラメーター設定手段と、
前記パラメーター設定手段により設定されたパラメーターを用いて、前記動態画像に対して２軸以上の特徴量空間でそれぞれの軸の特徴量でフィルタリングを施すフィルタリング手段と、
前記フィルタリング手段においてフィルタリングが施された動態画像に基づいて、前記検査対象部位に関する特徴量を算出する特徴量算出手段と、
を備え、
前記パラメーター設定手段は、前記２軸以上の特徴量空間のグラフを表示するグラフ表示手段を備え、前記グラフ表示手段により表示されたグラフ上におけるユーザー操作に基づいて、前記パラメーターを設定する動態解析装置。
【請求項２】
前記フィルタリング手段は、互いに異なる特性を持つ２軸以上の特徴量空間でフィルタリングを施す請求項１に記載の動態解析装置。
【請求項３】
前記フィルタリングに用いられる特徴量は、少なくとも一の前記動態画像の画素信号値の時間変化に関する特徴量を含む請求項１又は２に記載の動態解析装置。
【請求項４】
前記動態画像の画素信号値の時間変化に関する特徴量は、前記動態画像の画素信号値の時間周波数、前記動態画像の画素信号値の時間方向のプロファイル形状に関する特徴量の何れか一以上である請求項３に記載の動態解析装置。
【請求項５】
前記画素信号値の時間方向のプロファイル形状に関する特徴量は、前記プロファイル形状に基づいて算出されるピーク間積分、ピーク間距離、歪度、速度比、加速度比、又は前記プロファイル形状を機械学習の識別器に入力したときの出力結果である請求項４に記載の動態解析装置。
【請求項６】
前記フィルタリングに用いられる特徴量は、少なくとも一の前記動態画像の画素信号値の時間変化に関する特徴量及び少なくとも一の前記動態画像の画素信号値の空間的変化に関する特徴量を含む請求項１〜５の何れか一項に記載の動態解析装置。
【請求項７】
前記動態画像の画素信号値の空間的変化に関する特徴量は、前記動態画像の空間周波数、前記動態画像の各フレーム画像における画素信号値の分布を示す曲面の形状を表す曲率特徴量の何れか一以上である請求項６に記載の動態解析装置。
【請求項８】
前記グラフ表示手段は、前記動態画像に含まれる信号成分を前記２軸以上の特徴量空間のグラフ上で分類したガイド情報を併せて表示する請求項１〜７の何れか一項に記載の動態解析装置。
【請求項９】
前記特徴量算出手段は、前記フィルタリング手段によりフィルタリングが施された前記動態画像の各フレーム画像を複数の小領域に分割し、複数のフレーム画像の対応する小領域間における下記の（１）〜（７）の何れか一以上を前記検査対象部位の動態に関する特徴量として算出する請求項１〜８の何れか一項に記載の動態解析装置。
（１）時間的に隣接する２つのフレーム画像間の画素信号値の差分
（２）最大画素信号値と最小画素信号値の差分
（３）各フレーム画像の画素信号値の最小画素信号値からの差分
（４）各フレーム画像の画素信号値の最大画素信号値からの差分
（５）画素信号値の時間方向のプロファイル形状に基づいて算出されるピーク間積分
（６）画素信号値の時間方向のプロファイル形状に基づいて算出されるピーク間距離
（７）画素信号値の時間方向のプロファイル形状に基づいて算出される速度比
（８）画素信号値の時間方向のプロファイル形状に基づいて算出される加速度比
【請求項１０】
前記特徴量算出手段による算出結果を表示する結果表示手段を備える請求項１〜９の何れか一項に記載の動態解析装置。
【請求項１１】
前記動態画像に含まれる特定の構造物による画像信号成分を減弱する減弱手段を備える請求項１〜１０の何れか一項に記載の動態解析装置。
【請求項１２】
検査対象部位に放射線を照射して動態撮影することにより動態画像を取得する撮影手段と、
請求項１〜１１の何れか一項に記載の動態解析装置と、
を備える動態解析システム。
【請求項１３】
検査対象部位に放射線を照射して動態撮影することにより得られた動態画像に対して２軸以上の特徴量空間でそれぞれの軸の特徴量でフィルタリングを施すためのパラメーターを設定するパラメーター設定工程と、
前記パラメーター設定工程において設定されたパラメーターを用いて、前記動態画像に対して２軸以上の特徴量空間でそれぞれの軸の特徴量でフィルタリングを施すフィルタリング工程と、
前記フィルタリング工程においてフィルタリングが施された動態画像に基づいて、前記検査対象部位に関する特徴量を算出する特徴量算出工程と、
を含み、
前記パラメーター設定工程は、前記２軸以上の特徴量空間のグラフを表示手段に表示させ、前記表示手段に表示されたグラフ上におけるユーザー操作に基づいて、前記パラメーターを設定する動態解析方法。
【請求項１４】
コンピューターを、
検査対象部位に放射線を照射して動態撮影することにより得られた動態画像に対して２軸以上の特徴量空間でそれぞれの軸の特徴量でフィルタリングを施すためのパラメーターを設定するパラメーター設定手段、
前記パラメーター設定手段により設定されたパラメーターを用いて、前記動態画像に対して２軸以上の特徴量空間でそれぞれの軸の特徴量でフィルタリングを施すフィルタリング手段、
前記フィルタリング手段においてフィルタリングが施された動態画像に基づいて、前記検査対象部位に関する特徴量を算出する特徴量算出手段、
として機能させ、
前記パラメーター設定手段は、前記２軸以上の特徴量空間のグラフを表示手段に表示させ、前記表示手段に表示されたグラフ上におけるユーザー操作に基づいて、前記パラメーターを設定するプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、動態解析装置、動態解析システム、動態解析方法及びプログラムに関する。
【背景技術】
【０００２】
従来のフィルム／スクリーンや輝尽性蛍光体プレートを用いた放射線（Ｘ線）の静止画撮影及び診断に対し、ＦＰＤ（flat panel detector）等の半導体イメージセンサーを利用して検査対象部位の動態画像を撮影し、診断に応用する試みがなされるようになってきている。具体的には、半導体イメージセンサーの画像データの読取・消去の応答性の早さを利用し、半導体イメージセンサーの読取・消去のタイミングと合わせて放射源からパルス状の放射線を連続的に照射し、１秒間に複数回の撮影を行って、検査対象部位の動態を撮影する。撮影により取得された一連の複数枚の画像を順次表示することにより、医師は検査対象部位の一連の動きを観察することが可能となる。
【０００３】
ところで、肺の診断においては、肺の機能（換気機能や肺血流機能）が低下している箇所がないかを観察することは重要である。スパイロメーターによる換気機能測定や、血液ガス分析による肺血流機能測定は比較的簡便な手法として重用されているが、これらからは局所情報を知りえず、肺のどの部分が機能低下しているかまでは把握できない。一方で局所解析のためには、肺換気シンチグラムによる換気機能測定、肺血流シンチグラムによる肺血流機能測定を用いることができるが、放射性物質の吸引や血管注入が必要であり、被検者への負担が大きいという問題がある。
【０００４】
そこで、上記半導体イメージセンサーを用いた動態撮影と、その撮影データ解析により、上記検査を代替する技術が開示されている。
【０００５】
例えば、特許文献１には、息止め時の胸部動態画像における肺内の画素値が心拍動性により変化する性質を利用し、ピクセル値の変化情報を肺血流の情報とみなして、肺塞栓症や心臓疾患等の診断のために有効利用する技術が開示されている。
【０００６】
また、特許文献２には、呼吸状態下で撮影された胸部動態画像を構成する各フレーム画像の肺野領域を複数の領域に分割し、各領域内の画素信号値の時間変化を時間方向のハイパスフィルタでフィルタリングして各領域内の肺血流に関する特徴量を生成する技術が開示されている。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特許第５０９３７２７号公報
【特許文献２】特開２０１４−１２８６８７号公報
【発明の概要】
【発明が解決しようとする課題】
【０００８】
動態撮影において、健常者であれば息止め撮影は比較的容易に可能であるが、救急外来患者・手術患者・小児・老人などへの撮影など、息止め撮影困難な状況が医療現場には散見される。よって、特許文献１に記載のように、息止め時のみ精度良好な解析手法は医療現場においては実用的ではない。
【０００９】
一方、呼吸状態下の動態画像には解析不要信号（ノイズ）が息止め時より多く混在するため、それを精度良く分離する必要がある。例えば、特許文献２においては、画素信号値の時間変化をカットオフ周波数が固定値の時間方向のハイパスフィルタでフィルタリングすることにより、血流信号成分とそれ以外の信号成分を分離している。しかし、単体の特徴量の固定値をパラメーターとしてフィルタリングを行うと、被検者の個体差等の種々の要因により精度良くノイズが除去できない場合があった。
【００１０】
本発明の課題は、動態画像におけるノイズ除去精度を向上させることである。
【課題を解決するための手段】
【００１１】
上記課題を解決するため、請求項１に記載の発明の動態解析装置は、
検査対象部位に放射線を照射して動態撮影することにより得られた動態画像に対して２軸以上の特徴量空間でそれぞれの軸の特徴量でフィルタリングを施すためのパラメーターを設定するパラメーター設定手段と、
前記パラメーター設定手段により設定されたパラメーターを用いて、前記動態画像に対して２軸以上の特徴量空間でそれぞれの軸の特徴量でフィルタリングを施すフィルタリング手段と、
前記フィルタリング手段においてフィルタリングが施された動態画像に基づいて、前記検査対象部位に関する特徴量を算出する特徴量算出手段と、
を備え、
前記パラメーター設定手段は、前記２軸以上の特徴量空間のグラフを表示するグラフ表示手段を備え、前記グラフ表示手段により表示されたグラフ上におけるユーザー操作に基づいて、前記パラメーターを設定する。
【００１２】
請求項２に記載の発明は、請求項１に記載の発明において、
前記フィルタリング手段は、互いに異なる特性を持つ２軸以上の特徴量空間でフィルタリングを施す。
【００１３】
請求項３に記載の発明は、請求項１又は２に記載の発明において、
前記フィルタリングに用いられる特徴量は、少なくとも一の前記動態画像の画素信号値の時間変化に関する特徴量を含む。
【００１４】
請求項４に記載の発明は、請求項３に記載の発明において、
前記動態画像の画素信号値の時間変化に関する特徴量は、前記動態画像の画素信号値の時間周波数、前記動態画像の画素信号値の時間方向のプロファイル形状に関する特徴量の何れか一以上である。
【００１５】
請求項５に記載の発明は、請求項４に記載の発明において、
前記画素信号値の時間方向のプロファイル形状に関する特徴量は、前記プロファイル形状に基づいて算出されるピーク間積分、ピーク間距離、歪度、速度比、加速度比、又は前記プロファイル形状を機械学習の識別器に入力したときの出力結果である。
【００１６】
請求項６に記載の発明は、請求項１〜５の何れか一項に記載の発明において、
前記フィルタリングに用いられる特徴量は、少なくとも一の前記動態画像の画素信号値の時間変化に関する特徴量及び少なくとも一の前記動態画像の画素信号値の空間的変化に関する特徴量を含む。
【００１７】
請求項７に記載の発明は、請求項６に記載の発明において、
前記動態画像の画素信号値の空間的変化に関する特徴量は、前記動態画像の空間周波数、前記動態画像の各フレーム画像における画素信号値の分布を示す曲面の形状を表す曲率特徴量の何れか一以上である。
【００２３】
請求項８に記載の発明は、請求項１〜７の何れか一項に記載の発明において、
前記グラフ表示手段は、前記動態画像に含まれる信号成分を前記２軸以上の特徴量空間のグラフ上で分類したガイド情報を併せて表示する。
【００２８】
請求項９に記載の発明は、請求項１〜８の何れか一項に記載の発明において、
前記特徴量算出手段は、前記フィルタリング手段によりフィルタリングが施された前記動態画像の各フレーム画像を複数の小領域に分割し、複数のフレーム画像の対応する小領域間における下記の（１）〜（７）の何れか一以上を前記検査対象部位の動態に関する特徴量として算出する。
（１）時間的に隣接する２つのフレーム画像間の画素信号値の差分
（２）最大画素信号値と最小画素信号値の差分
（３）各フレーム画像の画素信号値の最小画素信号値からの差分
（４）各フレーム画像の画素信号値の最大画素信号値からの差分
（５）画素信号値の時間方向のプロファイル形状に基づいて算出されるピーク間積分
（６）画素信号値の時間方向のプロファイル形状に基づいて算出されるピーク間距離
（７）画素信号値の時間方向のプロファイル形状に基づいて算出される速度比
（８）画素信号値の時間方向のプロファイル形状に基づいて算出される加速度比
【００２９】
請求項１０に記載の発明は、請求項１〜９の何れか一項に記載の発明において、
前記特徴量算出手段による算出結果を表示する結果表示手段を備える。
【００３０】
請求項１１に記載の発明は、請求項１〜１０の何れか一項に記載の発明において、
前記動態画像に含まれる特定の構造物による画像信号成分を減弱する減弱手段を備える。
【００３１】
請求項１２に記載の発明の動態解析システムは、
検査対象部位に放射線を照射して動態撮影することにより動態画像を取得する撮影手段と、
請求項１〜１１の何れか一項に記載の動態解析装置と、
を備える。
【００３２】
請求項１３に記載の発明の動態解析方法は、
検査対象部位に放射線を照射して動態撮影することにより得られた動態画像に対して２軸以上の特徴量空間でそれぞれの軸の特徴量でフィルタリングを施すためのパラメーターを設定するパラメーター設定工程と、
前記パラメーター設定工程において設定されたパラメーターを用いて、前記動態画像に対して２軸以上の特徴量空間でそれぞれの軸の特徴量でフィルタリングを施すフィルタリング工程と、
前記フィルタリング工程においてフィルタリングが施された動態画像に基づいて、前記検査対象部位に関する特徴量を算出する特徴量算出工程と、
を含み、
前記パラメーター設定工程は、前記２軸以上の特徴量空間のグラフを表示手段に表示させ、前記表示手段に表示されたグラフ上におけるユーザー操作に基づいて、前記パラメーターを設定する。
【００３３】
請求項１４に記載の発明のプログラムは、
コンピューターを、
検査対象部位に放射線を照射して動態撮影することにより得られた動態画像に対して２軸以上の特徴量空間でそれぞれの軸の特徴量でフィルタリングを施すためのパラメーターを設定するパラメーター設定手段、
前記パラメーター設定手段により設定されたパラメーターを用いて、前記動態画像に対して２軸以上の特徴量空間でそれぞれの軸の特徴量でフィルタリングを施すフィルタリング手段、
前記フィルタリング手段においてフィルタリングが施された動態画像に基づいて、前記検査対象部位に関する特徴量を算出する特徴量算出手段、
として機能させ、
前記パラメーター設定手段は、前記２軸以上の特徴量空間のグラフを表示手段に表示させ、前記表示手段に表示されたグラフ上におけるユーザー操作に基づいて、前記パラメーターを設定する。
【発明の効果】
【００３４】
本発明によれば、動態画像におけるノイズ除去精度を向上させることが可能となる。
【図面の簡単な説明】
【００３５】
【図１】本発明の実施形態における動態解析システムの全体構成を示す図である。
【図２】図１の撮影用コンソールの制御部により実行される撮影制御処理を示すフローチャートである。
【図３】図１の診断用コンソールの制御部により実行される画像解析処理を示すフローチャートである。
【図４】横軸を動態画像の画素信号値の時間周波数、縦軸を動態画像の画素信号値の空間周波数とした２軸の特徴量空間のグラフ上において、胸部の動態画像に含まれる信号成分を分類した結果を示す図である。
【図５】、肺血流に伴う画素信号値の時間変化とノイズに伴う画素信号値の時間変化を示す図である。
【図６】ピーク間積分を説明するための図である。
【図７】速度比を説明するための図である。
【図８】加速度比を説明するための図である。
【図９】フィルタリングパラメーターを設定するための入力ウィンドウの一例を示す図である。
【図１０】フィルタリングパラメーターを設定するための入力ウィンドウの一例を示す図である。
【図１１】検査対象部位の動態画像に含まれ得る信号成分を２軸以上の特徴量空間のグラフ上に分類したガイド情報が表示されたフィルタリングパラメーターの入力ウィンドウの一例を示す図である。
【図１２】被検者背景情報と各特徴量の推奨パラメーターとを対応付けた表を表示したフィルタリングパラメーターの入力ウィンドウの一例を示す図である。
【図１３】撮影条件と各特徴量のパラメーターとを対応付けた表を表示したフィルタリングパラメーターの入力ウィンドウの一例を示す図である。
【図１４】心壁の動作周波数の算出手法を説明するための図である。
【図１５】横隔膜の動作周波数及び胸郭の動作周波数を算出する際の注目領域の設定を説明するための図である。
【図１６】縦隔主軸の動作周波数の算出手法を説明するための図である。
【発明を実施するための形態】
【００３６】
以下、図面を参照して本発明の実施形態を詳細に説明する。ただし、発明の範囲は、図示例に限定されない。
【００３７】
〔動態解析システム１００の構成 〕
まず、構成を説明する。
図１に、本実施形態における動態解析システム１００の全体構成を示す。
図１に示すように、動態解析システム１００は、撮影装置１と、撮影用コンソール２とが通信ケーブル等により接続され、撮影用コンソール２と、診断用コンソール３とがＬＡＮ（Local Area Network）等の通信ネットワークＮＴを介して接続されて構成されている。なお、本実施形態においては、撮影用コンソール２及び診断用コンソール３が通信ネットワークＮＴに有線により接続されている場合を例にとり説明するが、無線により接続されていることとしてもよい。
【００３８】
〔撮影装置１の構成〕
撮影装置１は、例えば、呼吸運動に伴う肺の膨張及び収縮の形態変化、心臓の拍動等の、周期性（サイクル）を持つ人体の動態を撮影する撮影手段である。動態撮影とは、被写体Ｍに対し、Ｘ線等の放射線をパルス状にして所定時間間隔で繰り返し照射するか(パルス照射)、もしくは、低線量率にして途切れなく継続して照射する(連続照射)ことで、動態を示す複数の画像を取得することをいう。動態撮影により得られた一連の画像を動態画像と呼ぶ。また、動態画像を構成する複数の画像のそれぞれをフレーム画像と呼ぶ。なお、以下の実施形態では、パルス照射により動態撮影を行う場合を例にとり説明する。
【００３９】
放射線源１１は、被写体Ｍを挟んで放射線検出部１３と対向する位置に配置され、放射線照射制御装置１２の制御に従って、被写体Ｍに対し放射線（Ｘ線）を照射する。
放射線照射制御装置１２は、撮影用コンソール２に接続されており、撮影用コンソール２から入力された放射線照射条件に基づいて放射線源１１を制御して放射線撮影を行う。撮影用コンソール２から入力される放射線照射条件は、例えば、パルスレート、パルス幅、パルス間隔、１撮影あたりの撮影フレーム数、Ｘ線管電流の値、Ｘ線管電圧の値、グリッドの種類、付加フィルタ種等である。パルスレートは、１秒あたりの放射線照射回数であり、後述するフレームレートと一致している。パルス幅は、放射線照射１回当たりの放射線照射時間である。パルス間隔は、１回の放射線照射開始から次の放射線照射開始までの時間であり、後述するフレーム間隔と一致している。
【００４０】
放射線検出部１３は、ＦＰＤ等の半導体イメージセンサーにより構成される。ＦＰＤは、例えば、ガラス基板等を有しており、基板上の所定位置に、放射線源１１から照射されて少なくとも被写体Ｍを透過した放射線をその強度に応じて検出し、検出した放射線を電気信号に変換して蓄積する複数の検出素子（画素）がマトリックス状に配列されている。各画素は、例えばＴＦＴ（Thin Film Transistor）等のスイッチング部を備えて構成されている。ＦＰＤにはＸ線をシンチレーターを介して光電変換素子により電気信号に変換する間接変換型、Ｘ線を直接的に電気信号に変換する直接変換型があるが、何れを用いてもよい。放射線検出部１３は、被写体Ｍを挟んで放射線源１１と対向するように設けられている。
【００４１】
ここで、放射線検出部１３において取得される画像の各画素の信号値（画素信号値）は、放射線検出部１３に到達した放射線強度を電気信号に変換した値、即ち、放射線検出部１３に到達した放射線強度に相関する値であり、到達した放射線強度が高いほど画素信号値は大きくなる。一方、本実施形態において、撮影用コンソール２及び診断用コンソール３では、各画素信号値を吸収線量を表すものとして取り扱う。即ち、撮影用コンソール２は、通信部２５により受信した一連のフレーム画像の各画素信号値を吸収線量を表す値に変換し、診断用コンソール３に送信する。即ち、撮影用コンソール２及び診断用コンソールでは、画素信号値は、吸収線量が多いほど画素信号値が高くなり、画像上では白く（低濃度で）描画される。
【００４２】
読取制御装置１４は、撮影用コンソール２に接続されている。読取制御装置１４は、撮影用コンソール２から入力された画像読取条件に基づいて放射線検出部１３の各画素のスイッチング部を制御して、当該各画素に蓄積された電気信号の読み取りをスイッチングしていき、放射線検出部１３に蓄積された電気信号を読み取ることにより、画像データを取得する。この画像データがフレーム画像である。そして、読取制御装置１４は、取得したフレーム画像を撮影用コンソール２に出力する。画像読取条件は、例えば、フレームレート、フレーム間隔、画素サイズ、画像サイズ（マトリックスサイズ）等である。フレームレートは、１秒あたりに取得するフレーム画像数であり、パルスレートと一致している。フレーム間隔は、１回のフレーム画像の取得動作開始から次のフレーム画像の取得動作開始までの時間であり、パルス間隔と一致している。
【００４３】
ここで、放射線照射制御装置１２と読取制御装置１４は互いに接続され、互いに同期信号をやりとりして放射線照射動作と画像の読み取りの動作を同調させるようになっている。
【００４４】
生体センサー１５は、動態撮影と同期した状態で、検査対象部位の動態に係る生体情報を取得して撮影用コンソール２に出力する。生体センサー１５は、例えば、パルスオキシメーター等を適用することができる。
【００４５】
〔撮影用コンソール２の構成〕
撮影用コンソール２は、放射線照射条件や画像読取条件を撮影装置１に出力して撮影装置１による放射線撮影及び放射線画像の読み取り動作を制御するとともに、撮影装置１により取得された動態画像を撮影技師等の撮影実施者によるポジショニングの確認や診断に適した画像であるか否かの確認用に表示する。
撮影用コンソール２は、図１に示すように、制御部２１、記憶部２２、操作部２３、表示部２４、通信部２５を備えて構成され、各部はバス２６により接続されている。
【００４６】
制御部２１は、ＣＰＵ（Central Processing Unit）、ＲＡＭ（Random Access Memory
）等により構成される。制御部２１のＣＰＵは、操作部２３の操作に応じて、記憶部２２に記憶されているシステムプログラムや各種処理プログラムを読み出してＲＡＭ内に展開し、展開されたプログラムに従って後述する撮影制御処理を始めとする各種処理を実行し、撮影用コンソール２各部の動作や、撮影装置１の放射線照射動作及び読み取り動作を集中制御する。
【００４７】
記憶部２２は、不揮発性の半導体メモリーやハードディスク等により構成される。記憶部２２は、制御部２１で実行される各種プログラムやプログラムにより処理の実行に必要なパラメーター、或いは処理結果等のデータを記憶する。例えば、記憶部２２は、図２に示す撮影制御処理を実行するためのプログラムを記憶している。また、記憶部２２は、検査対象部位（ここでは胸部）及び解析対象の種類に対応付けて放射線照射条件及び画像読取条件を記憶している。各種プログラムは、読取可能なプログラムコードの形態で格納され、制御部２１は、当該プログラムコードに従った動作を逐次実行する。
【００４８】
操作部２３は、カーソルキー、数字入力キー、及び各種機能キー等を備えたキーボードと、マウス等のポインティングデバイスを備えて構成され、キーボードに対するキー操作やマウス操作により入力された指示信号を制御部２１に出力する。また、操作部２３は、表示部２４の表示画面にタッチパネルを備えても良く、この場合、タッチパネルを介して入力された指示信号を制御部２１に出力する。
【００４９】
表示部２４は、ＬＣＤ（Liquid Crystal Display）やＣＲＴ（Cathode Ray Tube）等のモニターにより構成され、制御部２１から入力される表示信号の指示に従って、操作部２３からの入力指示やデータ等を表示する。
【００５０】
通信部２５は、ＬＡＮアダプターやモデムやＴＡ（Terminal Adapter）等を備え、通信ネットワークＮＴに接続された各装置との間のデータ送受信を制御する。
【００５１】
〔診断用コンソール３の構成〕
診断用コンソール３は、撮影用コンソール２から動態画像を取得し、取得した動態画像を解析し、その解析結果を表示して医師の診断を支援するための動態解析装置である。
診断用コンソール３は、図１に示すように、制御部３１、記憶部３２、操作部３３、表示部３４、通信部３５を備えて構成され、各部はバス３６により接続されている。
【００５２】
制御部３１は、ＣＰＵ、ＲＡＭ等により構成される。制御部３１のＣＰＵは、操作部３３の操作に応じて、記憶部３２に記憶されているシステムプログラムや、各種処理プログラムを読み出してＲＡＭ内に展開し、展開されたプログラムに従って、後述する画像解析処理を始めとする各種処理を実行し、診断用コンソール３各部の動作を集中制御する。制御部３１は、パラメーター設定手段、フィルタリング手段、特徴量算出手段、減弱手段として機能する。
【００５３】
記憶部３２は、不揮発性の半導体メモリーやハードディスク等により構成される。記憶部３２は、制御部３１で画像解析処理を実行するためのプログラムを始めとする各種プログラムやプログラムにより処理の実行に必要なパラメーター、或いは処理結果等のデータを記憶する。これらの各種プログラムは、読取可能なプログラムコードの形態で格納され、制御部３１は、当該プログラムコードに従った動作を逐次実行する。
【００５４】
また、記憶部３２は、後述する画像解析処理においてノイズ除去のためのフィルタリングに用いられたパラメーター（フィルタリングパラメーターと呼ぶ）の値や、性別や年齢等の被検者背景情報に応じた推奨パラメーターの値を記憶する。
【００５５】
操作部３３は、カーソルキー、数字入力キー、及び各種機能キー等を備えたキーボードと、マウス等のポインティングデバイスを備えて構成され、キーボードに対するキー操作やマウス操作により入力された指示信号を制御部３１に出力する。また、操作部３３は、表示部３４の表示画面にタッチパネルを備えても良く、この場合、タッチパネルを介して入力された指示信号を制御部３１に出力する。
【００５６】
表示部３４は、ＬＣＤやＣＲＴ等のモニターにより構成され、制御部３１から入力される表示信号の指示に従って、各種表示を行う。表示部３４は、グラフ表示手段、結果表示手段として機能する。
【００５７】
通信部３５は、ＬＡＮアダプターやモデムやＴＡ等を備え、通信ネットワークＮＴに接続された各装置との間のデータ送受信を制御する。
【００５８】
〔動態解析システム１００の動作〕
次に、上記動態解析システム１００における動作について説明する。
【００５９】
（撮影装置１、撮影用コンソール２の動作）
まず、撮影装置１、撮影用コンソール２による撮影動作について説明する。
図２に、撮影用コンソール２の制御部２１において実行される撮影制御処理を示す。撮影制御処理は、制御部２１と記憶部２２に記憶されているプログラムとの協働により実行される。
【００６０】
まず、撮影実施者により撮影用コンソール２の操作部２３が操作され、被写体Ｍとなる被検者の情報（被検者情報呼ぶ。例えば、被検者ＩＤ、氏名、被検者背景情報（身長、体重、年齢、性別等の被検者の身体的特徴等））、検査情報（検査対象部位（ここでは、胸部）、解析対象の種類（例えば、換気、肺血流）等）の入力が行われる（ステップＳ１）。なお、被検者情報や検査情報は、通信部２５を介してＲＩＳ（Radiology Information System）から受信する構成としてもよい。
【００６１】
次いで、入力された検査情報に基づいて、放射線照射条件が記憶部２２から読み出されて放射線照射制御装置１２に設定されるとともに、画像読取条件が記憶部２２から読み出されて読取制御装置１４に設定される（ステップＳ２）。なお、放射線照射条件及び画像読取条件は、操作部２３により撮影実施者が調整し、設定し直すことができる。
【００６２】
次いで、操作部２３の操作による放射線照射の指示が待機される（ステップＳ３）。ここで、例えば、呼吸状態下の撮影であれば、被検者（被写体Ｍ）に楽にするように指示し、安静呼吸を促す。撮影準備が整った時点で、操作部２３を操作して放射線照射指示を入力する。
【００６３】
操作部２３により放射線照射指示が入力されると（ステップＳ３；ＹＥＳ）、放射線照射制御装置１２及び読取制御装置１４に撮影開始指示が出力され、動態撮影が開始される（ステップＳ４）。即ち、放射線照射制御装置１２に設定されたパルス間隔で放射線源１１により放射線が照射され、放射線検出部１３によりフレーム画像が取得される。
【００６４】
予め定められたフレーム数の撮影が終了すると、制御部２１により放射線照射制御装置１２及び読取制御装置１４に撮影終了の指示が出力され、撮影動作が停止される。撮影されるフレーム数は、少なくとも１呼吸サイクルが撮影できる枚数である。
【００６５】
撮影により取得されたフレーム画像は順次撮影用コンソール２に入力され、撮影順を示す番号（フレーム番号）と対応付けて記憶部２２に記憶されるとともに（ステップＳ５）、表示部２４に表示される（ステップＳ６）。撮影実施者は、表示された動態画像によりポジショニング等を確認し、撮影により診断に適した画像が取得された（撮影ＯＫ）か、再撮影が必要（撮影ＮＧ）か、を判断する。そして、操作部２３を操作して、判断結果を入力する。
【００６６】
操作部２３の所定の操作により撮影ＯＫを示す判断結果が入力されると（ステップＳ７；ＹＥＳ）、動態撮影で取得された一連のフレーム画像のそれぞれに、動態画像を識別するための識別ＩＤや、被検者情報、検査情報、放射線照射条件、画像読取条件、撮影順を示す番号（フレーム番号）等の情報が付帯され（例えば、画像データのヘッダ領域に書き込まれ）、通信部２５を介して診断用コンソール３に送信される（ステップＳ８）。そして、本処理は終了する。一方、操作部２３の所定の操作により撮影ＮＧを示す判断結果が入力されると（ステップＳ７；ＮＯ）、記憶部２２に記憶された一連のフレーム画像が削除され（ステップＳ９）、本処理は終了する。この場合、再撮影が必要となる。
【００６７】
（診断用コンソール３の動作）
次に、診断用コンソール３における動作について説明する。
診断用コンソール３においては、通信部３５を介して撮影用コンソール２から動態画像の一連のフレーム画像が受信されると、制御部３１と記憶部３２に記憶されているプログラムとの協働により図３に示す画像解析処理が実行される。
【００６８】
以下、図３を参照して画像解析処理の流れについて説明する。
まず、ノイズ除去のためのフィルタリングに用いられる特徴量が設定される（ステップＳ１１）。
【００６９】
ここで、本願発明者らの検討によれば、動態画像に含まれる信号成分は、少なくとも互いに異なる特性を持つ２軸以上の特徴量空間、例えば、動態画像の画素信号値の時間変化に関する特徴量及び動態画像の画素信号値の空間的変化に関する特徴量を含む２軸以上の特徴量空間上で分類することができる。異なる特性を持つとは、同一特徴量ではないことを含む、互いに従属していないこと、を意味する。従属しているとは、例えばｙ＝３ｘのように、片方の変数が決定することでもう一方の変数の値が決まることを表す。同じ軸(例えば、空間周波数が２本)や、表しているものが本質的に同じ軸（例えば、半径と面積）は除外される。
【００７０】
図４に、横軸を動態画像の画素信号値の時間変化の周波数（以下、時間周波数と呼ぶ）、縦軸を動態画像の画素信号値の空間的変化の周波数（以下、空間周波数と呼ぶ）とした２軸の特徴量空間のグラフ上において、胸部の動態画像に含まれる信号成分を分類した結果を示す。即ち、互いに異なる特性を持つ２軸以上の特徴量空間のそれぞれの特徴量に、抽出したい信号成分の上限及び／下限をフィルタリングパラメーターとして設定してフィルタリングを行うことにより、それ以外の信号成分（即ちノイズ）を精度良く除去することが可能となる。
【００７１】
動態画像の画素信号値の時間変化に関する特徴量としては、例えば、上述の時間周波数（横軸を動態撮影開始からの経過時間、縦軸を画素信号値としてプロットしたときの波形の周波数）の他、動態画像の画素信号値の時間方向のプロファイル形状（上記の波形の形状）に関する特徴量が挙げられる。心拍動は、収縮期と拡張期では収縮期の方が相対的に短時間であり、肺血流に伴う画素信号値の時間変化は、図５に示すように、“急坂を登り、なだらかな坂を下りる”ように、非対称のピークとなる。これに対し、ノイズ成分の画素信号値の時間変化は、図５に示すように、そのような挙動を示さない。そこで、肺血流の信号成分（血流信号成分）を抽出したい場合は、動態画像の画素信号値の時間方向のプロファイル形状のピークの非対称性を表す特徴量をフィルタリングに用いる特徴量とすることができる。
【００７２】
動態画像の画素信号値の時間方向のプロファイル形状のピークの非対称性を表す特徴量としては、例えば、ピーク間積分、ピーク間距離、速度比、加速度比、歪度、上述のプロファイル形状を機械学習の識別器に入力したときに識別器から出力される血流らしさを示す値（尤度）等が挙げられる。
【００７３】
ピーク間積分は、各周期における信号値上昇時の信号積分値Supと、信号値下降時の信号積分値Sdown（図６参照）に関する統計値である。例えば、以下の（式１）〜（式３）の何れかにより求めることができる。ここで、jは、周期番号、MAXは、最大値、MEDは、中央値を示す（（式４）〜（式１２）についても同様）。
Σ Sdownj  / Σ Supj                   （式１）
MAX（Sdownj）/ MAX(Supj）              （式２）
MED（Sdownj）/ MED(Supj）              （式３）
【００７４】
ピーク間距離は、各周期における信号値上昇時の上昇時間Bupと、信号値下降時の下降時間Bdown（図６参照）に関する統計値である。例えば、以下の（式４）〜（式６）の何れかにより求めることができる。
Σ Bdownj / Σ Bupj                     (式４)
MAX (Bdownj）/ MAX（Bupj）             （式５）
MED (Bdownj）/ MED（Bupj）             （式６）
【００７５】
速度比は、各周期における傾きの最大値Vmaxと最小値Vmin（図７参照）に関する統計値である。例えば、以下の（式７）〜（式９）の何れかにより求めることができる。ABSは、絶対値を示す。
Σ ABS(Vmaxj) / Σ ABS(Vminj)          （式７）
MAX( ABS(Vmaxj)） / MAX( ABS(Vminj)）  （式８）
MED( ABS(Vmaxj)） / MED( ABS(Vminj)）  （式９）
上記の速度比、ピーク間積分、ピーク間差分は、血流信号では高値となり、ノイズは低値となる。
【００７６】
加速度比は、図８に示すように、各周期における画素信号値を微分し、その微分値につき、正値の領域(図８のドットで示すエリア)における最大微分値(信号値の２次微分)Amaxと、負値の領域(図８の斜線で示すエリア)における最小微分値(信号値の２次微分)Aminに関する以下の（式１０）〜（式１２）の何れかにより求めることができる。
Σ ABS(Amaxj) / Σ ABS(Aminj)          （式１０）
MAX( ABS(Amaxj) / MAX(ABS(Aminj)）     （式１１）
MED( ABS(Amaxj) / MED(ABS(Aminj)）     （式１２）
【００７７】
上述のピーク間積分、ピーク間差分、速度比、加速度比は、血流信号は高値となり、ノイズは低値となる。
【００７８】
歪度は、まず、各周期の歪度を以下の（式１３）により算出し、算出した各周期の歪度の統計値を、動態画像の画素信号値の時間変化に関する特徴量とすることができる。各周期の歪度の統計値は、例えば、各周期の歪度の最大値、上位から所望のパーセンテージの平均値、最頻値等である。例えば、最大値の場合、血流信号は高値となり、ノイズは低値となる。
【数１】
0006701880.tif 000002
【００７９】
また、動態画像の画素信号値の時間方向のプロファイル形状に関する特徴量としては、動態画像の画素信号値の時間方向のプロファイル形状を示すデータを機械学習の識別器に入力したときに識別器から出力される出力結果を用いることもできる。例えば、公知文献１、２に記載のディープラーニング、ＤＮＮ（Deep Natural Network）等に上記のプロファイル形状を示すデータを入力したときに出力される血流クラスの尤度等を用いることができる（公知文献１：岡谷貴之「ディープラーニングと画像認識―基礎と最近の動向―」，オペレーションズ・リサーチ 60(4), 191-197,2015、公知文献２：人工知能学会監修「深層学習」近代科学社(2015)）。例えば、動態画像における特定時間内の注目領域（例えば、５画素×５画素）の画素信号値を所定回数（例えば、２５回）サンプリングした元データ（例えば、２５×２５のデータ）を公知文献１に記載の畳込みニューラルネットワーク(Convolutional natural network)に入力したときに畳込みニューラルネットワークから出力される血流クラスの尤度（血流らしさを示す値）を用いることができる。
【００８０】
動態画像の画素信号値の空間的変化に関する特徴量としては、例えば、上述の空間周波数の他、曲率特徴量としてもよい。曲率特徴量は、曲線や曲面の曲がり具合を表す指標であり、画像の位置座標を示すx,y軸に対して画素信号値をｚ軸にプロットした３次元空間における形状（即ち、画像上の画素信号値の分布を示す曲面の形状）を数値表現できる。入力画像の各座標x,yに対応する画素信号値f(x,y)に対して曲率特徴量を計算し、その曲率特徴量に対してフィルタリングパラメーターを設定してフィルタリングを行うことで、上記の形状を絞り込むことが可能となる。
曲率特徴量としては、例えば、特開２００６−２３０９０４号公報に記載のShapeIndex値SI(x,y)を用いることができる。また、特開２００５−８０７５８号公報に記載のガウス曲率や平均曲率を用いることとしてもよい。例えば、曲率特徴量としてShapeIndex値SI(x,y)を用いる場合、SI(x,y)のフィルタリングパラメーターを０．２５未満に自動設定し、下記の数式によりフィルタリングを実施することができる。
・SI(x,y)＜0.25のとき、f(x,y)＝f(x,y)
・SI(x,y)≧0.25のとき、f(x,y)＝０
【００８１】
ステップＳ１１においては、上述の複数の特徴量等の中から、互いに異なる特性を持つ２軸以上の特徴量が設定される。特徴量の設定は、操作部３３によりユーザーが設定することができる。なお、フィルタリングに用いられる特徴量が予め決まっている（例えば、特徴量として時間周波数と空間周波数のみが用意されている）場合には、ステップＳ１１は省略することができる。
以下の説明では、時間周波数と空間周波数がフィルタリングに用いられる特徴量として設定された場合を例にとり説明する。
【００８２】
次いで、設定された特徴量のそれぞれにフィルタリングパラメーターが設定される（ステップＳ１２）。ステップＳ１２においては、ステップＳ１１で設定された特徴量のそれぞれについて、抽出範囲の上限及び／又は下限がフィルタリングパラメーターとして設定される。フィルタリングパラメーターの設定は、手動により行う構成としてもよいし、自動により行う構成としてもよい。
【００８３】
以下、ステップＳ１２におけるフィルタリングパラメーターの設定の例について説明する。
例えば、図９に示すような各特徴量の上限値及び下限値の入力欄が設けられた入力ウィンドウ３４１が表示部３４に表示され、ユーザーによる操作部３３の操作により入力された数値がフィルタリングパラメーターとして設定される。
【００８４】
または、図１０に示すように、設定された特徴量のそれぞれを軸とした特徴量空間のグラフが表示された入力ウィンドウ３４２が表示部３４に表示され、ユーザーによる操作部３３の操作により各軸の上限位置及び下限位置が指定されると、指定された位置に対応する値がフィルタリングパラメーターとして設定される。例えば、特徴量空間のグラフに各軸の上限位置を示す線、下限位置を示す線が表示され、各線を操作部３３によりスライドさせることにより、各軸の上限位置及び下限位置を指定することができる。また、各軸には、「上限開放」「下限開放」を指定するためのチェックボックスが表示され、チェックボックスをチェックすることで、上限値や下限値を開放することができる。この場合、チェックされた値に対応する線は非表示となる。例えば、図１０に示すように、時間周波数の上限開放にチェックが行われた場合は、時間周波数の上限位置を示す線が非表示となる。
また、図１１に示すように、入力ウィンドウ３４２において、検査対象部位の動態画像に含まれ得る信号成分を２軸以上の特徴量空間のグラフ上で分類したガイド情報Ｇを併せて表示することとしてもよい。これにより、ユーザーは特徴量の持つ意味合いを直感的に理解することができ、容易にフィルタリングパラメーターの値を設定することが可能となる。
なお、図９〜図１１では、フィルタリングパラメーターの設定範囲が１箇所である場合を例示したが、２箇所以上に設定可能な構成としてもよい。
【００８５】
または、記憶部３２に記憶されている、過去検査において設定されたフィルタリングパラメーターを参照することにより、フィルタリングパラメーターが設定される。
例えば、同一被検者の過去検査において設定されたフィルタリングパラメーター（今回の動態画像に付帯されている付帯情報と同じ被検者ＩＤ、同じ解析対象の種類に対応付けて記憶部３２に記憶されている、ステップＳ１１で設定された各特徴量に対するフィルタリングパラメーター）が記憶部３２から検索され、検索されたフィルタリングパラメーターの値が今回のフィルタリングパラメーターとして設定される。
【００８６】
または、性別・年齢等の被検者背景情報と各特徴量の推奨するフィルタリングパラメーター（推奨パラメーターと呼ぶ）とを対応付けた表（対応表）を記憶部３２に記憶しておき、動態画像に付帯されている付帯情報の被検者背景情報に対応付けられている推奨パラメーターのうち、解析対象の種類及び設定された特徴量に応じた推奨パラメーターの値が今回のフィルタリングパラメーターとして設定される。この場合、例えば、図１２に示すように、記憶部３２に記憶されている対応表を表示部３４に表示し、ユーザーによる操作部３３の操作に応じてフィルタリングパラメーターを設定することとしてもよい。
【００８７】
または、撮影条件（放射線照射条件）と過去にその撮影条件で撮影された動態画像に対して設定された各特徴量のフィルタリングパラメーターの実績（過去実績）とを対応付けた表（対応表）を記憶部３２に記憶しておき、記憶部３２に記憶されている対応表において、動態画像に付帯されている付帯情報の撮影条件と同じ条件に対応付けられている、設定された特徴量のフィルタリングパラメーターの値が今回のフィルタリングパラメーターとして設定される。この場合、例えば、図１３に示すように、記憶部３２に記憶されている対応表を表示部３４に表示し、ユーザーによる操作部３３の操作によって選択されたフィルタリングパラメーターを設定することとしてもよい。また、図１３に示すように、過去にそのフィルタリングパラメーターが選択された（設定された）回数を併せて表示することとしてもよい。これにより、ユーザーが容易にフィルタリングパラメーターを設定することが可能となる。
また、撮影条件と過去実績とを対応付けた対応表を被検者背景情報（例えば、年齢、性別等）毎に記憶部３２に記憶しておき、記憶部３２に記憶されている、動態画像に付帯されている付帯情報の被検者背景情報に対応する対応表において、動態画像に付帯されている付帯情報の撮影条件と同じ条件に対応付けられている、設定された特徴量のフィルタリングパラメーターの値が今回のフィルタリングパラメーターとして設定されるようにしてもよい。或いは、記憶部３２に記憶されている、付帯情報の被検者背景情報（年齢、性別等）に対応する対応表を表示部３４に表示し、ユーザーによる操作部３３の操作によって選択されたフィルタリングパラメーターを設定することとしてもよい。即ち、表示部３４に表示するフィルタリングパラメーター候補を被検者背景情報により絞り込めるようにしてもよい。例えば、４０代の女性が被検者である場合、４０代の女性に限った対応表を表示部３４に表示することとしてもよい。これにより、被検者の背景が略同一の集団での過去実績を鑑みて撮影条件を選択することができるため、より精度良くフィルタリングパラメーターを設定することができる。
【００８８】
または、動態画像に含まれる構造物の動態を表す特徴量が算出され、算出された特徴量に基づいて、自動的に、動態画像の画素信号値の時間変化に関する特徴量（ここでは、時間周波数）についてのフィルタリングパラメーターが設定される。
例えば、心壁の動作及び肺血流は、心臓の拍動によって生じるものであるため、心壁の動作周波数は、動態画像における肺血流による信号成分の時間周波数と略一致しているとみなすことができる。また、横隔膜の移動や、胸郭や体表皮膚境界面の移動は、肺の呼吸運動（換気）により生じるものであるため、横隔膜、胸郭、体表皮膚境界面の動作周波数は、動態画像における換気による信号成分の時間周波数と略一致しているとみなすことができる。また、縦隔主軸は、体動があると動くため、縦隔主軸の動作周波数は、動態画像における体動による信号成分の時間周波数と略一致しているとみなすことができる。
そこで、ステップＳ１２においては、例えば、動態画像を解析して心壁、横隔膜（又は、胸郭、体表皮膚境界面）、縦隔主軸のそれぞれの動作周期λを求め、求めた動作周期λに基づいて動作周波数ｆを求め、求めた動作周波数ｆに基づいて、動態画像の画素信号値の時間周波数のフィルタリングパラメーターが設定される。
【００８９】
ここで、心壁の動作周波数ｆは、以下のようにして求めることができる。まず、フレーム画像毎に、心臓輪郭の抽出を行う。心臓輪郭の抽出は、公知の画像処理技術、例えば、特許第２７９６３８１号公報に記載の心臓輪郭決定方法等を用いて行うことができる。次いで、図１４に示すように、左心室により形成される左第４弓に基準点Ｐを設定する。例えば、抽出された心臓輪郭の高さＨを３等分し、下から１／３の高さの水平線を測定水平線ｌとし、その測定水平線ｌと抽出された左心室の心臓輪郭との交点を基準点Ｐとする。そして、基準点Ｐのｘ座標を時系列にプロットすることにより得られる波形の隣り合う頂点の時刻差から周期λを求め、周波数ｆ＝１／λより、心壁の動作周波数ｆを算出する。なお、肺血流観測のための時間周波数把握の観点からは本来は肺へ血液を送り込む右心室の壁を観察することが望ましいが、それと左心室の壁とのそれぞれの移動周期は充分近似すると考え、安定的に自動取得できる左心室の壁の輪郭を計算に用いている。
また、抽出された心臓輪郭の横断面の直径Ｄを時系列にプロットすることにより得られる波形の隣り合う頂点の時刻差から周期λを求め、周波数ｆ＝１／λより、心壁の動作周波数ｆを求めることとしてもよい。
【００９０】
横隔膜の動作周波数ｆは、以下のようにして求めることができる。まず、フレーム画像毎に、横隔膜境界抽出を行う。横隔膜境界抽出は、公知の画像処理技術、例えば、特表２００２−５０３８６１号公報に記載の反復全体しきい値分析と局所しきい値分析をベースとした手法等を用いて行うことができる。横隔膜境界抽出後、図１５（ａ）に示すように、横隔膜境界部のある垂直方向断面の位置に領域Ｒ１を設定し、この領域Ｒ１における横隔膜境界のｙ座標の代表値（例えば、平均値、中央値等）を時系列にプロットすることにより得られる波形の隣り合う頂点の時刻差から周期λを求める。または、横隔膜境界部に図１５（ｂ）に示すように領域Ｒ１１を設定し、この領域Ｒ１１内の画素信号値の代表値（例えば、平均値、中央値）を時系列にプロットすることにより得られる波形の隣り合う頂点の時刻差から周期λを求めることとしてもよい。周波数ｆ＝１／λより、横隔膜の動作周波数ｆを算出することができる。なお、画素信号値を利用する場合は、設定する領域Ｒ１１は横隔膜境界付近であれば全フレーム画像において必ずしも横隔膜境界が入っていなくてもよい。
【００９１】
胸郭の動作周波数ｆは、以下のようにして求めることができる。まず、フレーム画像毎に、胸郭境界抽出を行う。胸郭境界抽出は、公知の画像処理技術、例えば、上記特表２００２−５０３８６１号公報に記載の手法等を用いて行うことができる。胸郭境界抽出後、図１５（ａ）に示すように、胸郭境界部のある水平方向断面の位置に領域Ｒ２を設定し、この領域Ｒ２における胸郭境界部のｘ座標の代表値（例えば、平均値、中央値等）を時系列にプロットすることにより得られる波形の隣り合う頂点の時刻差から周期λを求める。または、胸郭境界部に図１５（ｂ）に示すように領域Ｒ１２を設定し、この領域Ｒ１２内の画素信号値の代表値（例えば、平均値、中央値）を時系列にプロットすることにより得られる波形の隣り合う頂点の時刻差から周期λを求めることとしてもよい。周波数ｆ＝１／λより、胸郭の動作周波数ｆを算出することができる。
【００９２】
体表皮膚境界面の動作周波数ｆは、以下のようにして求めることができる。まず、フレーム画像毎に、体表皮膚境界面の抽出を行う。体表皮膚境界面の抽出は、公知の画像処理技術、例えば、２値化処理等を用いて行うことができる。体表皮膚境界面の抽出後、図１５（ａ）に示すように、体表皮膚境界面のある水平方向断面の位置に領域Ｒ３を設定し、この領域Ｒ３における体表皮膚境界面のｘ座標の代表値（例えば、平均値、中央値等）を時系列にプロットすることにより得られる波形の隣り合う頂点の時刻差から周期λを求める。または、体表皮膚境界面に図１５（ｂ）に示すように領域Ｒ１３を設定し、この領域Ｒ１３内の画素信号値の代表値（例えば、平均値、中央値）を時系列にプロットすることにより得られる波形の隣り合う頂点の時刻差から周期λを求めることとしてもよい。周波数ｆ＝１／λより、体表皮膚境界面の動作周波数ｆを算出することができる。
【００９３】
縦隔主軸の動作周波数ｆは、以下のようにして求めることができる。まず、フレーム画像毎に、肺野領域の抽出を行う。肺野領域の抽出は、公知の画像処理技術、例えば、上記特表２００２−５０３８６１号公報に記載の手法等を用いて行うことができる。肺野領域の抽出後、図１６に示すように、まず、右肺野領域及び左肺野領域の夫々の輪郭線において、ｘ軸方向においては、全輪郭線のy座標最小値(最も頭部に近い側)の１点から正中線までの間に位置する範囲、y軸方向においては、全輪郭線のy軸方向の長さＬ（y座標の最大値-最小値）の所望の割合(例えば１／３)分の長さを算出し、y座標最小値から算出した長さまでの範囲を計算対象となる線分（図１６の白色線）を定める。次いで、計算対象となった線分において、その線分の重心の位置を算出し、当該重心の位置の時間プロファイル（時間変化を示すグラフ）を作成する。そして、重心がフレーム画像のある基準位置から指定距離以上離れ、かつ、指定距離以下まで近づく軌跡をもつ現象を検出した場合、それにかかる時間を体動による揺らぎの周期λとし、周波数ｆ=１／λより、縦隔主軸の動作周波数ｆ（体動の周波数ｆ）を算出することができる。
また、体動の代表ケースは主に左右の揺らぎであることから、前記重心に代わり、計算対象となった線分の各点のｘ座標の代表値(平均値、中央値、最頻値、など)を求めて時間プロファイル（時間変化を示すグラフ）を作成し、振幅が所定以上の波が検出された場合、その波のピーク間の時間を体動による揺らぎの周期λとし、周波数ｆ=１／λより、縦隔主軸の動作周波数ｆ（体動の周波数ｆ）を算出することができる。
【００９４】
例えば、解析対象が肺血流である場合、心壁の動作周波数ｆを含み、横隔膜（又は、胸郭、体表皮膚境界面）及び縦隔主軸の動作周波数ｆを含まない範囲が抽出されるように時間周波数のフィルタリングパラメーターが設定される。例えば、解析対象が換気である場合、横隔膜（又は、胸郭、体表皮膚境界面）の動作周波数ｆを含み、心壁及び縦隔主軸の動作周波数ｆを含まない範囲が抽出されるように時間周波数のフィルタリングパラメーターが設定される。
【００９５】
または、生体センサー１５により検出された生体情報に基づいて、自動的に動態画像の画素信号値の時間変化に関する特徴量についてのフィルタリングパラメーターを設定することとしてもよい。例えば、生体センサー１５がパルスオキシメーターの場合、生体情報として、動脈血酸素飽和度（ＳｐＯ２）や脈拍数（ＰＲ）を取得することができる。また、吸気時には胸腔内圧が下がり中心静脈圧が下がることで静脈還流が増加し、これによりパルスオキシメーターが装着されている末梢の静脈床の容積がわずかに減少するため、パルスオキシメーターにより検出される脈波のベースラインが下がる。呼気時にはその反対の現象が生じる。これを元に、呼吸数を推定することができる（Addison PS, et al. Developing an algorithm for pulse oximetry derived respiratory rate (RRoxi): a healthy volunteer study. J Clin Monit Comput. 2012;26(1):45-51.）。そこで、生体センサー１５から得られる脈拍数、呼吸数のそれぞれの逆数を求めることで、脈拍や呼吸の周波数を求め、求めた脈拍や呼吸の周波数に基づいて、動態画像の画素信号値の時間変化に関する特徴量についてのフィルタリングパラメーターを設定することができる。
【００９６】
上記説明した何れの手法によりフィルタリングパラメーターを設定するかは、操作部３３によりユーザーが設定することができる。
【００９７】
フィルタリングパラメーターの設定が終了すると、各フレーム画像にローカルマッチング処理及びワーピング処理（非線形歪変換処理）が施され、フレーム画像間で肺野領域の位置合わせが行われる（ステップＳ１３）。例えば、各フレーム画像において、同等の解剖学的特徴位置の複数の点（ランドマーク）が抽出され、１番目に撮影されたフレーム画像を基準画像として、基準画像からのランドマーク各点のシフト値が算出される。そして、算出されたランドマークの各点のシフト値に基づいて、全画素のシフト値がｎ次元多項式の近似処理を用いて算出され、算出されたシフト値分だけ各画素をシフトした画像が生成される。
なお、本実施形態においては、安静呼吸時に撮影を行っているが、安静呼吸時においては呼吸運動による肺野領域の位置ずれはわずかである。そのため、ステップＳ１３は省略することとしてもよい。
【００９８】
次いで、構造物の減弱処理が行われる（ステップＳ１４）。ステップＳ１４においては、各フレーム画像から解析に不要な所定の構造物の画像信号成分の減弱処理が行われる。以下、所定の構造物が肋骨である場合について説明する。
まず、各フレーム画像から肋骨の抽出が行われる。肋骨の抽出は、公知の画像処理技術を適用することができる。例えば、予め用意した肋骨テンプレートとのテンプレートマッチングや、エッジ検出、エッジ検出後にカーブフィッティング関数を当てはめる等の手法を適用することができる。また、特開２００５−２０３３８号公報に提案されているように、エッジ検出（放物線近似形状検出）により検出された肋骨の初期形状を肋骨形状モデル（教師データから得られた平均形状と、教師データの主成分分析により得られた複数の主成分形状との線形和として任意の肋骨形状を生成する）に投影して、肋骨のモデル投影形状を求める手法を適用することとしてもよい。
各フレーム画像から肋骨が抽出されると、各フレーム画像毎に、抽出された各肋骨領域断面の画素値プロファイル（横軸をフレーム画像における垂直方向の位置、縦軸を画素信号値としてプロットしたもの）が作成され、作成されたプロファイルにローパスフィルタを適用することによりノイズ等の空間的な高周波成分が除去された後、ノイズ等が除去された画素値プロファイルの値がもとのフレーム画像から減算される。これにより、肋骨による画像信号成分を減弱したフレーム画像を得ることができる。
なお、本実施形態においては、各フレーム画像から肋骨を減弱する場合を例にとり説明したが、例えば、心臓、縦隔、横隔膜等を減弱することとしてもよい。また、ステップＳ１５のフィルタリング処理後に構造物の減弱処理を行うこととしてもよい。
【００９９】
次いで、設定された特徴量空間（ここでは、２軸の特徴量空間）の設定されたフィルタリングパラメーターを用いて、動態画像にノイズを除去するためのフィルタリング処理が行われる（ステップＳ１５）。例えば、設定された特徴量が動態画像の画素信号値の時間変化に関する特徴量及び空間的変化に関する特徴量である場合、まず、動態画像の画素信号値の空間的変化に関する特徴量によるフィルタリング処理が実行され、次いで、動態画像の画素信号値の時間変化に関する特徴量によるフィルタリング処理が実行される。例えば、設定された特徴量が時間周波数と空間周波数である場合、まず、フレーム画像毎に、フーリエ変換が行われ、設定された空間周波数のパラメーターに基づくバンドパスフィルタ、ハイパスフィルタ、またはローパスフィルタでフィルタリングが行われ、逆フーリエ変換される。或いは、各フレーム画像の肺野領域を、設定されたフィルタリングパラメーターの空間周波数に応じたサイズの複数の小領域に分割し、各小領域内において、画素信号値の代表値（例えば、平均値、中央値等）を算出し、領域内の画素信号値を算出した代表値に置き換えることで空間周波数によるフィルタリングを行うこととしてもよい。次いで、空間周波数でフィルタリング済みの一連のフレーム画像の対応する画素（又は、小領域）の画素信号値が時間方向にプロットされることにより取得される画素信号値の時間変化が、設定された時間周波数のフィルタリングパラメーターに応じたカットオフ周波数の時間方向のハイパスフィルタ、ローパスフィルタ、またはバンドパスフィルタでフィルタリングされる。これにより、動態画像におけるノイズ除去を精度良く行うことができる。
【０１００】
次いで、フィルタリングに用いられた各特徴量のパラメーターが記憶部３２に記憶される（ステップＳ１６）。例えば、フィルタリングに用いられた各特徴量のパラメーターが今回の動態画像の付帯情報に対応付けて記憶部３２に記憶される。また、撮影条件（放射線照射条件）と過去に設定されたフィルタリングパラメーターの実績（過去実績）とを対応付けた表（対応表）が記憶部３２に記憶されている場合は、対応表から今回使用された撮影条件及び今回使用された各特徴量のフィルタリングパラメーターが一致するレコードが検索され、一致するレコードが検索された場合は、選択回数が１インクリメントされる。一致するレコードが検索されなかった場合は、対応表に新たなレコードが追加され、今回使用された撮影条件及び今回使用された各特徴量のフィルタリングパラメーターの組み合わせが格納され、選択回数に１が設定される。
【０１０１】
次いで、フィルタリング処理後の一連のフレーム画像の肺野領域が解析のための複数の小領域に分割され、一連のフレーム画像の対応する小領域毎に、肺野の動態に関する特徴量が算出される（ステップＳ１７）。例えば、小領域毎に、画素信号値の代表値（例えば、平均値、中央値等）の時間変化が算出され、この時間変化に基づいて、以下の（１）〜（８）の何れか一以上の特徴量が算出される。なお、解析のための小領域は、１画素単位としてもよい。
（１）時間的に隣接する２つのフレーム画像間の画素信号値の差分（フレーム間差分）
（２）最大画素信号値と最小画素信号値の差分、
（３）各フレーム画像の画素信号値の最小画素信号値からの差分
（４）各フレーム画像の画素信号値の最大画素信号値からの差分
（５）画素信号値の時間方向のプロファイル形状に基づいて算出されるピーク間積分
（６）画素信号値の時間方向のプロファイル形状に基づいて算出されるピーク間距離
（７）画素信号値の時間方向のプロファイル形状に基づいて算出される速度比
（８）画素信号値の時間方向のプロファイル形状に基づいて算出される加速度比
【０１０２】
次いで、ステップＳ１７における特徴量の算出結果が表示部３４に表示され（ステップＳ１８）、画像解析処理は終了する。
例えば、ステップＳ１７において、フレーム間差分、各フレーム画像の画素信号値の最小画素信号値からの差分、各フレーム画像の画素信号値の最大画素信号値からの差分の何れかが算出された場合は、各フレーム画像の各小領域を特徴量の値に応じた輝度若しくは色で示して、一連のフレーム画像を表示部３４に動画表示又は並べて表示する。ステップＳ１７において最大画素信号値と最小画素信号値の差分、ピーク間積分、ピーク間距離、速度比、加速度比の何れかが算出された場合は、各小領域を特徴量の値に応じた輝度若しくは色で示して、表示部３４に表示することとしてもよい。または、ステップＳ１７において、上述の（式１）〜（式１２）の統計記号（Σ、MAX、MED）を抜いた式により周期毎に算出される値を特徴量として算出し、各周期毎の各小領域を特徴量に応じた輝度若しくは色で示した周期毎の画像を表示部３４に動画表示することとしてもよい。
【０１０３】
以上説明したように、診断用コンソール３によれば、制御部３１は、検査対象部位に放射線を照射して動態撮影することにより得られた動態画像に対して２軸以上の特徴量空間でそれぞれの軸の特徴量でフィルタリングを施し、フィルタリングが施された動態画像に基づいて、検査対象部位の動態に関する特徴量を算出する。従って、動態画像におけるノイズ除去精度を向上させることができる。
【０１０４】
例えば、制御部３１は、動態画像に対して２軸以上の特徴量空間のそれぞれの軸の特徴量でフィルタリングを施すためのパラメーターを設定し、設定されたパラメーターを用いて動態画像に対してフィルタリングを施す。
従って、動態画像毎に個別に複数軸の特徴量空間のフィルタリングパラメーターを設定するため、複数の特徴量に対する個体差に応じたパラメーター設定が可能となり、動態画像におけるノイズ除去精度を向上させることができる。特に、呼吸状態下の胸部動態画像にはノイズが息止め時よりも多く混在するため、その効果はより顕著である。また、動態画像から精度良くノイズ除去を行うことができるため、その後の解析が失敗してやり直しとなることを防止することができ、医療現場における作業効率を向上させることができる。また、解析の精度も向上させることができる。
【０１０５】
また、例えば、制御部３１は、フィルタリングに用いる２軸以上の特徴量空間のグラフを表示部３４に表示し、表示されたグラフ上におけるユーザー操作に基づいて、フィルタリングに用いるパラメーターを設定する。従って、ユーザーが直感的に理解しやすいＧＵＩで複数の特徴量のパラメーターを設定することができる。また、動態画像に含まれ得る信号成分を２軸以上の特徴量空間のグラフ上で分類したガイド情報を併せて表示することで、ユーザーが特徴量の持つ意味合いを直感的に理解することが可能となるので、より容易に、ユーザーが納得してパラメーターを設定することが可能となる。
【０１０６】
また、例えば、制御部３１は、同一被検者の過去検査において設定されたパラメーターを取得し、取得したパラメーターに基づいて、フィルタリングに用いるパラメーターを設定する。従って、ユーザーの手間を必要とせずにパラメーターを設定することが可能となる。
【０１０７】
また、例えば、制御部３１は、動態画像に付帯されている被検者背景情報に対応する推奨パラメーターを記憶部３２から取得し、取得した推奨パラメーターに基づいて、フィルタリングに用いるパラメーターを設定する。従って、被検者の性別や年齢等の特徴に応じたパラメーターを容易に設定することが可能となる。
【０１０８】
また、制御部３１は、動態画像と略一致する撮影条件で撮影された動態画像に対して過去に設定されたパラメーターを記憶部３２から取得し、取得したパラメーターに基づいて、フィルタリングに用いるパラメーターを設定する。従って、撮影された撮影条件に応じたパラメーターを容易に設定することが可能となる。
【０１０９】
また、制御部３１は、動態画像を解析して当該動態画像に含まれる所定の構造物の動態を表す特徴量を取得し、取得した所定の構造物の動態を表す特徴量に基づいて、動態画像の画素信号値の時間変化に関する特徴量についてのパラメーターを設定する。従って、動態画像に最適なパラメーターをユーザーの手間を必要とせずに自動的に設定することが可能となる。
【０１１０】
また、制御部３１は、生体センサー１５から取得された生体情報に基づいて、動態画像の画素信号値の時間変化に関する特徴量についてのパラメーターを設定する。従って、被検者の検査対象部位の動態に応じた最適なパラメーターをユーザーの手間を必要とせずに自動的に設定することが可能となる。
【０１１１】
また、制御部３１は、動態画像に含まれる特定の構造物による画像信号成分を減弱した後で動態画像にフィルタリングを施すので、予め構造物によるノイズを動態画像から除去しておくことができ、より一層精度良くノイズ除去を行うことが可能となる。
【０１１２】
また、制御部３１は、小領域毎に算出した肺野の動態に関する特徴量を表示部３４に表示するので、ユーザーは、小領域毎に、肺野の動態に関する特徴量を把握することが可能となる。
【０１１３】
なお、本実施形態における記述は、本発明に係る好適な動態解析システムの一例であり、これに限定されるものではない。
【０１１４】
例えば、上記実施形態においては、胸部を検査対象部位とした場合を例にとり説明したが、他の部位を検査対象部位としてもよい。
【０１１５】
また、例えば、上記の説明では、本発明に係るプログラムのコンピューター読み取り可能な媒体としてハードディスクや半導体の不揮発性メモリー等を使用した例を開示したが、この例に限定されない。その他のコンピューター読み取り可能な媒体として、ＣＤ-ＲＯＭ等の可搬型記録媒体を適用することが可能である。また、本発明に係るプログラムのデータを通信回線を介して提供する媒体として、キャリアウエーブ(搬送波)も適用される。
【０１１６】
その他、動態解析システム１００を構成する各装置の細部構成及び細部動作に関しても、本発明の趣旨を逸脱することのない範囲で適宜変更可能である。
【符号の説明】
【０１１７】
１００  動態解析システム
１  撮影装置
１１  放射線源
１２  放射線照射制御装置
１３  放射線検出部
１４  読取制御装置
２  撮影用コンソール
２１  制御部
２２  記憶部
２３  操作部
２４  表示部
２５  通信部
２６  バス
３  診断用コンソール
３１  制御部
３２  記憶部
３３  操作部
３４  表示部
３５  通信部
３６  バス
【図１】
0006701880.tif 000003
【図２】
0006701880.tif 000004
【図３】
0006701880.tif 000005
【図４】
0006701880.tif 000006
【図５】
0006701880.tif 000007
【図６】
0006701880.tif 000008
【図７】
0006701880.tif 000009
【図８】
0006701880.tif 000010
【図９】
0006701880.tif 000011
【図１０】
0006701880.tif 000012
【図１１】
0006701880.tif 000013
【図１２】
0006701880.tif 000014
【図１３】
0006701880.tif 000015
【図１４】
0006701880.tif 000016
【図１５】
0006701880.tif 000017
【図１６】
0006701880.tif 000018