知財求人 - 知財ポータルサイト「IP Force」
▶ ラドウィスプ プライベート リミテッドの特許一覧 ▶ 株式会社メディオットの特許一覧
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2022095871
(43)【公開日】2022-06-28
(54)【発明の名称】診断支援プログラム
(51)【国際特許分類】
A61B 6/00 20060101AFI20220621BHJP
【ＦＩ】
A61B6/00 330A
A61B6/00 360Z
A61B6/00 350S
【審査請求】有
【請求項の数】6
【出願形態】ＯＬ
(21)【出願番号】P 2022066013
(22)【出願日】2022-04-13
(62)【分割の表示】P 2019563986の分割
【原出願日】2019-01-07
(31)【優先権主張番号】P 2018000855
(32)【優先日】2018-01-05
(33)【優先権主張国・地域又は機関】JP
(31)【優先権主張番号】P 2018043511
(32)【優先日】2018-03-09
(33)【優先権主張国・地域又は機関】JP
(71)【出願人】
【識別番号】517441952
【氏名又は名称】パラメヴィア プライベート リミテッド
(71)【出願人】
【識別番号】516216335
【氏名又は名称】株式会社メディオット
(74)【代理人】
【識別番号】100114258
【弁理士】
【氏名又は名称】福地 武雄
(74)【代理人】
【識別番号】100125391
【弁理士】
【氏名又は名称】白川 洋一
(72)【発明者】
【氏名】阿部 武彦
(72)【発明者】
【氏名】吉田 典史
(57)【要約】 （修正有）
【課題】呼気または吸気の全部または一部を含む呼吸要素ごとに形状が変化する領域の動きを表示することが可能な診断支援プログラムを提供する。
【解決手段】各フレーム画像の特定領域の画素に基づいて、呼気または吸気の全部または一部を含む呼吸要素の周期を特定する処理と、特定した呼吸要素の周期に基づいて、肺野を検出する処理と、検出した肺野を複数のブロック領域に分割し、各フレーム画像におけるブロック領域の画像の変化を計算する処理と、各フレーム画像における各ブロック領域の画像の変化をフーリエ変換する処理と、フーリエ変換後に得られるスペクトルのうち、呼吸要素の周期に対応するスペクトルを含む一定の帯域内のスペクトルを抽出する処理と、前記一定の帯域から抽出したスペクトルに対し逆フーリエ変換する処理と、逆フーリエ変換後の各画像をディスプレイに表示する処理と、を含む。
【選択図】図６Ａ
【解決手段】各フレーム画像の特定領域の画素に基づいて、呼気または吸気の全部または一部を含む呼吸要素の周期を特定する処理と、特定した呼吸要素の周期に基づいて、肺野を検出する処理と、検出した肺野を複数のブロック領域に分割し、各フレーム画像におけるブロック領域の画像の変化を計算する処理と、各フレーム画像における各ブロック領域の画像の変化をフーリエ変換する処理と、フーリエ変換後に得られるスペクトルのうち、呼吸要素の周期に対応するスペクトルを含む一定の帯域内のスペクトルを抽出する処理と、前記一定の帯域から抽出したスペクトルに対し逆フーリエ変換する処理と、逆フーリエ変換後の各画像をディスプレイに表示する処理と、を含む。
【選択図】図６Ａ
【特許請求の範囲】
【請求項1】
人体の画像を解析し、解析結果を表示する診断支援プログラムであって、
複数のフレーム画像を取得する処理と、
前記各フレーム画像内の画素値を相対値または対数値として表示する処理と、
いずれかのフレーム画像と任意の間隔で選択された他のフレーム画像、または前記相対値若しくは対数値として表示されたいずれかのフレーム画像と、前記相対値若しくは対数値として表示され、任意の間隔で選択された他のフレーム画像との差分または比を算出する処理と、
前記算出された差分または比を可視化する処理と、
を備え、
肺野の濃度変化または血管の位置関係に基づいて、一つ以上のフレーム画像の肺野領域の画素値を算出することを特徴とする診断支援プログラム。
複数のフレーム画像を取得する処理と、
前記各フレーム画像内の画素値を相対値または対数値として表示する処理と、
いずれかのフレーム画像と任意の間隔で選択された他のフレーム画像、または前記相対値若しくは対数値として表示されたいずれかのフレーム画像と、前記相対値若しくは対数値として表示され、任意の間隔で選択された他のフレーム画像との差分または比を算出する処理と、
前記算出された差分または比を可視化する処理と、
を備え、
肺野の濃度変化または血管の位置関係に基づいて、一つ以上のフレーム画像の肺野領域の画素値を算出することを特徴とする診断支援プログラム。
【請求項2】
人体の画像を解析し、解析結果を表示する診断支援プログラムであって、
複数のフレーム画像を取得する処理と、
前記各フレーム画像内の画素値を相対値または対数値として表示する処理と、
いずれかのフレーム画像と任意の間隔で選択された他のフレーム画像、または前記相対値若しくは対数値として表示されたいずれかのフレーム画像と前記相対値若しくは対数値として表示され任意の間隔で選択された他のフレーム画像との差分または比を算出する処理と、
前記算出された差分または比を可視化する処理と、
を備え、
アーティファクトを除外したフレーム画像の画素値を算出することを特徴とする診断支援プログラム。
複数のフレーム画像を取得する処理と、
前記各フレーム画像内の画素値を相対値または対数値として表示する処理と、
いずれかのフレーム画像と任意の間隔で選択された他のフレーム画像、または前記相対値若しくは対数値として表示されたいずれかのフレーム画像と前記相対値若しくは対数値として表示され任意の間隔で選択された他のフレーム画像との差分または比を算出する処理と、
前記算出された差分または比を可視化する処理と、
を備え、
アーティファクトを除外したフレーム画像の画素値を算出することを特徴とする診断支援プログラム。
【請求項3】
心室容量の波形に基づいて、心臓または肺の動的な画像を表示することを特徴とする請求項１または請求項２記載の診断支援プログラム。
【請求項4】
肺野を複数のブロック領域に分割し、各フレーム画像におけるブロック領域の平均デンシティを算出することを特徴とする請求項１または請求項２に記載の診断支援プログラム。
【請求項5】
各フレーム画像のデンシティの変化を取得する際に、フィルタ処理を行なうことを特徴とする請求項１または請求項２に記載の診断支援プログラム。
【請求項6】
前記可視化する処理がされた信号をディスプレイに表示する処理をさらに有し、
前記信号をカラー表示し、または画素値の割合の値として表示することを特徴とする請求項１または請求項２に記載の診断支援プログラム。
前記信号をカラー表示し、または画素値の割合の値として表示することを特徴とする請求項１または請求項２に記載の診断支援プログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、人体の画像を解析し、解析結果を表示する技術に関する。
【背景技術】
【0002】
医師が胸部の動態画像によって肺の診断を行なう際、被写体が自然呼吸の状態で撮影された時系列の胸部動態画像の観察が重要である。生理学的なデータを取得しやすいスパイロメータ、ＲＩ（Radio Isotope）検査、形態的なデータが得られる単純Ｘ線写真、ＣＴ（Computed Tomography）などが肺機能を評価するための手法として知られている。しかし、生理学的なデータと形態的なデータの両者を効率良く取得することは容易ではない。
【0003】
近年、ＦＰＤ（Flat panel detector）等の半導体イメージセンサを利用して、人体の胸部の動態画像を撮影し、診断に用いる方法が試みられている。例えば、非特許文献１には、動態画像を構成する複数のフレーム画像の間で、信号値の差を示す差分画像を生成し、その差分画像から各信号値の最大値を求めて表示する技術が開示されている。
【0004】
また、特許文献１には、人体の胸部の動態を示す複数のフレーム画像のそれぞれのフレーム画像から肺野領域を抽出し、その肺野領域を複数の小領域に分割し、複数のフレーム画像の間において、分割された小領域を互いに対応付けて解析する技術が開示されている。この技術によれば、分割された小領域の動きを示す特徴量が表示される。
【先行技術文献】
【特許文献】
【0005】
【特許文献1】特許第５８７４６３６号明細書
【非特許文献】
【0006】
【非特許文献1】“Basic Imaging Properties of a Large Image Intensifier-TV Digital Chest Radiographic System” Investigative Radiology:April 1987; 22 : 328-335.
【発明の概要】
【発明が解決しようとする課題】
【0007】
しかしながら、非特許文献１に記載されている技術のように、単に動態画像の画素ごとのフレーム間差分値の最大値を表示するだけでは、医師が病態を把握することは容易ではない。また、特許文献１に記載されている技術のように、特徴量を表示するだけでは、やはり病態の把握には十分ではない。このため、呼吸や肺血管の状態に即した画像を表示することが望ましい。すなわち、被写体である人体の呼吸状態および血管動態全体を把握し、呼吸、心臓、肺門部の血管または血流の波形若しくは周波数、または画像の変化傾向に基づいて、実際の動きを示す画像を表示することが望ましい。
【0008】
本発明は、このような事情に鑑みてなされたものであり、呼気または吸気の全部または一部を含む呼吸要素ごとに形状が変化する領域の動きを表示することが可能な診断支援プログラムを提供することを目的とする。より具体的には、計測しようとしている新たな対象のデータに対し、既に取得している波の形およびHzに対する一致率やその他の不一致率を数値化することによって、診断の補助となる数値を算出し、さらに、これらの数値を画像化することにより、診断の補助となる画像を生成することを目的とする。
【課題を解決するための手段】
【0009】
（１）上記の目的を達成するために、本願は、以下のような手段を講じた。すなわち、本発明の一態様に係る診断支援プログラムは、人体の画像を解析し、解析結果を表示する診断支援プログラムであって、前記画像を格納するデータベースから複数のフレーム画像を取得する処理と、前記各フレーム画像の特定領域の画素に基づいて、呼気または吸気の全部または一部を含む呼吸要素の少なくとも一つの周波数を特定する処理と、前記特定した呼吸要素の少なくとも一つの周波数に基づいて、肺野を検出する処理と、前記検出した肺野を複数のブロック領域に分割し、前記各フレーム画像におけるブロック領域の画像の変化を計算する処理と、前記各フレーム画像における各ブロック領域の画像の変化をフーリエ変換する処理と、前記フーリエ変換後に得られるスペクトルのうち、前記呼吸要素の少なくとも一つの周波数に対応するスペクトルを含む一定の帯域内のスペクトルを抽出する処理と、前記一定の帯域から抽出したスペクトルに対して逆フーリエ変換する処理と、前記逆フーリエ変換後の各画像をディスプレイに表示する処理と、をコンピュータに実行させることを特徴とする。
【0010】
（２）また、本発明の一態様に係る診断支援プログラムは、前記フーリエ変換後に得られるスペクトルのうち、ノイズの周波数を含み、前記フレーム画像から得られる呼吸要素の周波数以外の周波数、または入力された周波数若しくは周波数帯域に対応するスペクトルを含む一定の帯域内のスペクトルを、フィルタを用いて抽出する処理をさらに含むことを特徴とする。
【0011】
（３）また、本発明の一態様に係る診断支援プログラムは、前記呼吸要素の周波数および前記各フレーム画像に基づいて、前記フレーム間の画像を生成する処理と、をさらに含むことを特徴とする。
【0012】
（４）また、本発明の一態様に係る診断支援プログラムは、人体の画像を解析し、解析結果を表示する診断支援プログラムであって、前記画像を格納するデータベースから複数のフレーム画像を取得する処理と、被写体の心拍または血管拍から抽出される心血管拍要素の少なくとも一つの周波数を特定する処理と、前記各フレーム画像の特定領域の画素に基づいて、呼気または吸気の全部または一部を含む呼吸要素の少なくとも一つの周波数を特定する処理と、前記特定した呼吸要素の少なくとも一つの周波数に基づいて、肺野を検出する処理と、前記検出した肺野を複数のブロック領域に分割し、前記各フレーム画像におけるブロック領域の画像の変化を計算する処理と、前記各フレーム画像における各ブロック領域の画像の変化をフーリエ変換する処理と、前記フーリエ変換後に得られるスペクトルのうち、前記心血管拍要素の少なくとも一つの周波数に対応するスペクトルを含む一定の帯域内のスペクトルを抽出する処理と、前記一定の帯域から抽出したスペクトルに対して逆フーリエ変換する処理と、前記逆フーリエ変換後の各画像をディスプレイに表示する処理と、をコンピュータに実行させることを特徴とする。
【0013】
（５）また、本発明の一態様に係る診断支援プログラムは、人体の画像を解析し、解析結果を表示する診断支援プログラムであって、前記画像を格納するデータベースから複数のフレーム画像を取得する処理と、被写体の心拍または血管拍から抽出される心血管拍要素の少なくとも一つの周波数を特定する処理と、肺野を検出する処理と、前記検出した肺野を複数のブロック領域に分割し、前記各フレーム画像におけるブロック領域の画像の変化を計算する処理と、前記各フレーム画像における各ブロック領域の画像の変化をフーリエ変換する処理と、前記フーリエ変換後に得られるスペクトルのうち、前記心血管拍要素の少なくとも一つの周波数に対応するスペクトルを含む一定の帯域内のスペクトルを抽出する処理と、前記一定の帯域から抽出したスペクトルに対して逆フーリエ変換する処理と、前記逆フーリエ変換後の各画像をディスプレイに表示する処理と、をコンピュータに実行させることを特徴とする。
【0014】
（６）また、本発明の一態様に係る診断支援プログラムは、前記フーリエ変換後に得られるスペクトルのうち、ノイズの周波数を含み、前記フレーム画像から得られる心血管拍要素の周波数以外の周波数、または入力された周波数若しくは周波数帯域に対応するスペクトルを含む一定の帯域内のスペクトルを、フィルタを用いて抽出する処理をさらに含むことを特徴とする。
【0015】
（７）また、本発明の一態様に係る診断支援プログラムは、前記特定した心血管拍要素の周波数および前記各フレーム画像に基づいて、前記フレーム間の画像を生成する処理をさらに含むことを特徴とする。
【0016】
（８）また、本発明の一態様に係る診断支援プログラムは、人体の画像を解析し、解析結果を表示する診断支援プログラムであって、前記画像を格納するデータベースから複数のフレーム画像を取得する処理と、被写体の血管拍から抽出される血管拍要素の少なくとも一つの周波数を特定する処理と、前記各フレーム画像について設定された解析範囲を複数のブロック領域に分割し、前記各フレーム画像におけるブロック領域の画像の変化を計算する処理と、前記各フレーム画像における各ブロック領域の画像の変化をフーリエ変換する処理と、前記フーリエ変換後に得られるスペクトルのうち、前記血管拍要素の少なくとも一つの周波数に対応するスペクトルを含む一定の帯域内のスペクトルを抽出する処理と、前記一定の帯域から抽出したスペクトルに対して逆フーリエ変換する処理と、前記逆フーリエ変換後の各画像をディスプレイに表示する処理と、をコンピュータに実行させることを特徴とする。
【0017】
（９）また、本発明の一態様に係る診断支援プログラムは、前記フーリエ変換後に得られるスペクトルのうち、ノイズの周波数を含み、前記フレーム画像から得られる血管拍要素の周波数以外の周波数、または入力された周波数若しくは周波数帯域に対応するスペクトルを含む一定の帯域内のスペクトルを、フィルタを用いて抽出する処理をさらに含むことを特徴とする。
【0018】
（１０）また、本発明の一態様に係る診断支援プログラムは、前記特定した心血管拍要素の周波数および前記各フレーム画像に基づいて、前記フレーム間の画像を生成する処理をさらに含むことを特徴とする。
【0019】
（１１）また、本発明の一態様に係る診断支援プログラムは、人体の画像を解析し、解析結果を表示する診断支援プログラムであって、前記画像を格納するデータベースから複数のフレーム画像を取得する処理と、前記各フレーム画像の特定領域の画素に基づいて、呼気または吸気の全部または一部を含む呼吸要素の少なくとも一つの周波数を特定する処理と、前記特定した呼吸要素の少なくとも一つの周波数に基づいて、肺野および横隔膜を検出する処理と、前記検出した肺野を複数のブロック領域に分割し、前記各フレーム画像におけるブロック領域の画素の変化率を算出する処理と、前記ブロック領域の画素の変化率と、呼吸と連動する動的部位の変化率との比の値である同調率を用いて、前記同調率が予め定められた一定の範囲内にあるブロック領域のみを抽出する処理と、前記抽出したブロック領域のみを含む各画像をディスプレイに表示する処理と、をコンピュータに実行させることを特徴とする。
【0020】
（１２）また、本発明の一態様に係る診断支援プログラムは、被写体の心拍または血管拍から抽出される心血管拍要素の少なくとも一つの周波数または血管拍から抽出される血管拍要素の少なくとも一つの周波数を特定する処理をさらに含むことを特徴とする。
【0021】
（１３）また、本発明の一態様に係る診断支援プログラムは、前記同調率の対数の値が、０を含む一定の範囲として定められることを特徴とする。
【0022】
（１４）また、本発明の一態様に係る診断支援プログラムは、特定のフレームにおいて検出した肺野上の少なくとも一つ以上のベジエ曲線（Bezier curve）を用いて、他のフレームにおける肺野を検出する処理をさらに含むことを特徴とする。
【0023】
（１５）また、本発明の一態様に係る診断支援プログラムは、前記検出した肺野内に内部制御点を選定し、前記肺野内の内部制御点を通る曲線または直線によって前記肺野を分割することを特徴とする。
【0024】
（１６）また、本発明の一態様に係る診断支援プログラムは、前記検出した肺野の外延およびその近傍における制御点の間隔を相対的に大きくし、前記検出した肺野内における部位毎の膨張比率に応じて、前記内部制御点の間隔を相対的に小さくすることを特徴とする。
【0025】
（１７）また、本発明の一態様に係る診断支援プログラムは、前記検出した肺野において、制御点の間隔を、人体に対して頭尾方向に進むに従って相対的に大きくし、または特定のベクトル方向に従って相対的に大きくすることを特徴とする。
【0026】
（１８）また、本発明の一態様に係る診断支援プログラムは、特定のフレームにおいて検出した肺野上の少なくとも一つ以上のベジエ曲面（Bezier surface）を用いて、他のフレームにおける肺野を検出する処理をさらに含むことを特徴とする。
【0027】
（１９）また、本発明の一態様に係る診断支援プログラムは、特定のフレームにおいて予め定められた解析範囲上に、少なくとも一つ以上のベジエ曲線（Bezier curve）を用いて、他のフレームにおいて前記解析範囲に対応する範囲を検出する処理をさらに含むことを特徴とする。
【0028】
（２０）また、本発明の一態様に係る診断支援プログラムは、少なくとも一つ以上のベジエ曲線（Bezier curve）を用いて、少なくとも肺野、血管または心臓を描画する処理をさらに含むことを特徴とする。
【0029】
（２１）また、本発明の一態様に係る診断支援プログラムは、人体の画像を解析し、解析結果を表示する診断支援プログラムであって、前記画像を格納するデータベースから複数のフレーム画像を取得する処理と、前記取得したすべてのフレーム画像について、ベジエ曲線を用いて解析範囲を特定する処理と、前記解析範囲内のインテンシティ（intensity）の変化に基づいて解析対象を検出する処理と、をコンピュータに実行させることを特徴とする。
【0030】
（２２）また、本発明の一態様に係る診断支援プログラムは、前記検出した解析対象の辺縁の特徴を算出する処理をさらに含むことを特徴とする。
【0031】
（２３）また、本発明の一態様に係る診断支援プログラムは、連続する各画像について、インテンシティ（intensity）の差分を算出することで横隔膜を検出し、前記検出した横隔膜または呼吸と連動する動的部位の位置または形状を示す指標を表示することを特徴とする。
【0032】
（２４）また、本発明の一態様に係る診断支援プログラムは、インテンシティ（intensity）の閾値を変化させることで、横隔膜以外の部位によって遮られていない横隔膜を表示し、横隔膜の全体形状を補間することを特徴とする。
【0033】
（２５）また、本発明の一態様に係る診断支援プログラムは、前記検出した横隔膜の位置若しくは形状、または呼吸と連動する動的部位の位置若しくは形状から、前記呼吸要素の少なくとも一つの周波数を計算する処理と、をさらに含むことを特徴とする。
【0034】
（２６）また、本発明の一態様に係る診断支援プログラムは、前記検出した肺野を空間的に正規化し、またはリコンストラクション（reconstruction）を利用して時間的に正規化する処理をさらに含むことを特徴とする。
【0035】
（２７）また、本発明の一態様に係る診断支援プログラムは、前記呼吸要素の少なくとも一つの周波数の位相を変化させ、または呼吸要素の波形を円滑化させることで、呼吸要素を補正することを特徴とする。
【0036】
（２８）また、本発明の一態様に係る診断支援プログラムは、解析範囲内のいずれかの部位の波形を特定し、前記特定した波形の周波数の構成要素を抽出し、前記波形の周波数の構成要素に対応する画像を出力することを特徴とする。
【0037】
（２９）また、本発明の一態様に係る診断支援プログラムは、解析範囲のデンシティ（density）を検出し、デンシティが相対的に大きく変化する箇所を除去することを特徴とする。
【0038】
（３０）また、本発明の一態様に係る診断支援プログラムは、前記フーリエ変換後に得られるスペクトルから、臓器特有の周期的な変化のスペクトル構成比に基づいて、逆フーリエ変換を行なう際の少なくとも一つの周波数を選択する処理をさらに含むことを特徴とする。
【0039】
（３１）また、本発明の一態様に係る診断支援プログラムは、前記呼吸要素の少なくとも一つの周波数に応じて、Ｘ線の照射間隔を調整するよう、Ｘ線撮影装置を制御することを特徴とする。
【0040】
（３２）また、本発明の一態様に係る診断支援プログラムは、前記逆フーリエ変換後に、振幅値が相対的に大きいブロックのみを抽出して表示することを特徴とする。
【0041】
（３３）また、本発明の一態様に係る診断支援プログラムは、前記肺野を同定した後、横隔膜または胸郭を特定し、横隔膜または胸郭の変化量を算出し、前記変化量から変化率を算出する処理をさらに含むことを特徴とする。
【0042】
（３４）また、本発明の一態様に係る診断支援プログラムは、特定のスペクトルに係数を乗算する処理をさらに含み、前記係数が乗算された特定にスペクトルに基づいて強調表示を行なうことを特徴とする。
【0043】
（３５）また、本発明の一態様に係る診断支援プログラムは、画像を格納するデータベースから複数のフレーム画像を取得した後、呼吸要素の周波数または波形を特定するために、解析対象となる部位にデジタルフィルタを施すことを特徴とする。
【0044】
（３６）また、本発明の一態様に係る診断支援プログラムは、前記各フレーム画像の特定領域の画素に基づいて、呼気または吸気の全部または一部を含む呼吸要素の複数の周波数を特定し、前記呼吸要素の複数の周波数のそれぞれに対応する各画像をディスプレイに表示することを特徴とする。
【0045】
（３７）また、本発明の一態様に係る診断支援プログラムは、ある一枚以上のフレーム画像の特定の範囲について、ある一定の値に集簇する画像を選択し、ディスプレイに表示することを特徴とする。
【発明の効果】
【0046】
本発明の一態様によれば、呼気または吸気の全部または一部を含む呼吸要素ごとに形状が変化する領域の動きを表示することが可能となる。
【図面の簡単な説明】
【0047】
【図1A】本実施形態に係る診断支援システムの概略構成を示す図である。
【図1B】肺領域の分割手法の一例を示す図である。
【図1C】肺の形態が時間の経過によって変化する様子を示す図である。
【図1D】肺の形態が時間の経過によって変化する様子を示す図である。
【図2A】特定ブロックの「intensity」変化と、それをフーリエ解析した結果を示す図である。
【図2B】心拍に近い周波数成分を抜き出したフーリエ変換結果と、これを逆フーリエ変換して心拍に近い周波数成分の「intensity」変化を示す図である。
【図2C】フーリエ変換後に得られたスペクトルのうち、ある一定の帯域を抽出する例を示す図である。
【図2D】肺の変化率を模式的に表した図である。
【図2E】肺野領域のパターン画像の例を示す図である。
【図2F】肺野領域のパターン画像の例を示す図である。
【図2G】肺野領域のパターン画像の例を示す図である。
【図2H】肺野領域のパターン画像の例を示す図である。
【図3A】肺野の輪郭を、ベジエ曲線および直線の両方を用いて描画した例を示す図であり、肺野が最大の状態を示す。
【図3B】肺野の輪郭を、ベジエ曲線および直線の両方を用いて描画した例を示す図であり、肺野が最小の状態を示す。
【図4A】前と次のフレームの間において、肺野の画像の前後を重ねた図である。
【図5】曲線回帰を行ない、横隔膜の相対的位置を近似した結果を示す図である。
【図6A】本実施形態に係る呼吸機能解析の概要を示すフローチャートである。
【図6B】ディスプレイに表示される画像の一例を示す図である。
【図6C】ディスプレイに表示される画像の一例を示す図である。
【図7】本実施形態に係る肺血流解析の概要を示すフローチャートである。
【図8】本実施形態に係るその他の血流解析の概要を示すフローチャートである。
【図9】フーリエ変換後に得られたスペクトルのうち、ある一定のスペクトルに係数を掛ける例を示す図である。
【図10】肺野を、ベジエ曲線を用いて描画した例を示す図である。
【図11】肺野を、ベジエ曲線を用いて分割した例を示す図である。
【図12】肺野を、ベジエ曲線を用いて分割した例を示す図である。
【図13】大動脈血流量の波形と心室容量の波形とを対比した一例を示す図である。
【図14】肺と肺の近辺の画素値の一例を示す図である。
【図15】人体の血管の概略構成を模式化した図である。
【発明を実施するための形態】
【0048】
まず、本発明の基本的な概念について説明する。本発明では、人体における呼吸や血管、肺野の面積及び体積、その他の生体運動において、一定の周期で反復するように捉えられる動きに対し、全体若しくはある部分的な範囲について、時間軸における一定の反復若しくは一定の運き（ルーチーン）を、波として捉え、計測する。波の計測結果については、（ア）波の形態自体、若しくは、（イ）波の間隔（周波数：Hz）を用いる。この２つの概念をまとめて「ベースデータ」と呼称する。
【0049】
同時期に同じようにリンクされる波が存在し得る。例えば、呼吸であれば、以下の近似を概念することができる。
（ある大雑把な範囲の「density」変化の平均）≒（胸郭の変化）≒（横隔膜の動き）≒（肺機能検査）≒（胸腹呼吸センサ）
（ある大雑把な範囲の「density」変化の平均）≒（胸郭の変化）≒（横隔膜の動き）≒（肺機能検査）≒（胸腹呼吸センサ）
【0050】
上記「（ア）波の形態自体」について、「波形同調性」という概念を用い、これに基づいて画像を表示する（Wave form tunable imaging）。また、上記「（イ）波の間隔（周波数：Hz）」について、「周波数同調性」という概念を用い、これに基づいて画像を表示する（Frequency tunable imaging）。
【0051】
例えば、心臓の場合、図１３に示す「大動脈血流量の波形と心室容量の波形とを対比した一例」のように、大動脈血流量のピークと、心室容積のピークや波形は一致しない。しかし、図１３中、時刻ｔ１からｔ２、時刻ｔ２からｔ３、時刻ｔ３からｔ４…といったように等間隔の時間的幅を１サイクルと定めると、大動脈血流量の１サイクルおよび心室容量の１サイクルが何度も繰り返されることとなり、各波形は周波数が同調していると言える。この波形について着目すると、図１３に示すような実測値から１サイクルを特定し、モデル波形を利用することによって、波の形（Wave form）を予想することができる。すなわち、「ベースデータとしての波形」の作り方として、実測しても良いし、周波数（サイクル）から生成しても良いし、モデル波形を利用しても良いし、個人間の波形を平均化して利用しても良い。心臓などの周波数を有する臓器のサイクル（周期）が分かれば、波の形（Wave form）が予想できるため、大動脈血流量や心室容量などの波形を把握し、この波形に基づいて、臓器の動的な画像を表示することが可能となる。
【0052】
なお、呼吸、心臓、肺門などの「density」の変化を取得する際、他の要素が混ざらないよう、予めデジタルフィルタをかけても良い。
【0053】
また、本発明では、「呼吸要素」という概念を用いる。「呼吸要素」とは、呼気または吸気の全部または一部を含む。例えば、「１呼吸」を、「１呼気」と「１吸気」に分けて考えることもできるし、「１呼気または１吸気」の「０％、１０％、２０％、３０％、４０％、５０％、６０％、７０％、８０％、９０％、１００％」のいずれかに限定して考えることもできる。さらに、各呼気の一定の割合のみ、例えば、呼気の１０％のみを抽出して評価することも可能である。これらのいずれかのデータ、またはこれらを組み合わせたデータを用いて、より精度の高い画像の抽出を可能とする。この際、何度も相互的に計算することもある。
【0054】
このような考え方は、「呼吸要素」のみならず、「心血管要素」についても同様に適用することが可能である。
【0055】
ここで、ベースデータを作る際、単一若しくは複数のモダリティから得られる特徴量（例えば、ある一定範囲の「density」、「volumetry」より構成される変化量、胸郭の動き、横隔膜の動き、「spirometry」、胸腹呼吸センサの2つ以上）、または、同じ呼吸周回などの複数回の波形測定により、お互いの成分抽出を補足し合い、精度を高めていく。これにより、アーティファクトの軽減、ライン（line）などのある一定の予想をもとに精度を高めることが可能となる。ここで、「density」とは、「密度」と訳されるが、画像においては、特定の領域における画素の「吸収値」を意味する。例えば、ＣＴでは、空気は「－１０００」、骨は「１０００」、水は「０」として用いられている。
【0056】
また、お互いの成分抽出による波の軸、幅、範囲およびHzの揺らぎ、幅を推定する。すなわち、複数回の重ね合わせによって、Hzの軸設定が平均化、分散によって軸、幅、範囲、Hzの最適レンジ（range）が計算される。その際に他の行動のHz（ノイズ）が抽出され、その波があればそれが入らない程度も相対的に計測していく場合がある。すなわち、波形要素の全体のうち、一部の波形のみを抽出する場合がある。
【0057】
本明細書では、「density」と「intensity」とを区別して用いる。「density」は、上述したように、吸収値を意味し、ＸＰやＸＰ動画の元画において、空気の透過性が高く、透過性の高い部分を白であることを数値にし、空気を「－１０００」、水を「０」、骨を「１０００」として表示するものとする。一方、「intensity」は、「density」から相対的に変化したもの、例えば、normalizedして濃度の幅、信号の程度に“変換”して表示したものとする。すなわち、「intensity」は、画像において、明暗や強調度などの相対的な値である。ＸＰ画像の吸収値を直接扱っている間は、「density」または「densityの変化（Δdensity）」として表す。そして、これを画像表現上の都合で、上記のような変換を行なって、「intensity」として表す。例えば、０から２５５の２５６階調にカラー表示する場合は、「intensity」となる。このような用語の区別は、ＸＰやＣＴの場合にあてはまる。
【0058】
一方、ＭＲＩの場合は、空気を「－１０００」、水を「０」、骨を「１０００」と定めようとしても、ＭＲＩの画素値、測定機械の種類、測定時の人の体調、体つき、測定時間によって、値が非常に変化してしまうという事情があり、また、Ｔ１強調像などＭＲＩの信号の採り方についても、その施設、測定機械の種類によって様々であり、一定ではない。このため、ＭＲＩの場合は、ＸＰやＣＴの場合のような「density」定義ができない。このため、ＭＲＩでは、最初に描出する段階から、相対値を取り扱うこととしており、最初から「intensity」として表現する。そして、その処理する信号も「intensity」である。
【0059】
以上により、ベースデータを得ることが可能となる。上記ベースデータに対し、計測したい新たな対象について、上記ベースデータの波形、波のHzのある一定の幅、範囲で抽出する。例えば、呼吸抽出のみや、血管抽出程度の幅、範囲、波形要素で抽出する。なお、この波形、Hzの幅に関しては、他の機能における波形要素、ノイズなどの「artifact」、他の同調性があると思われる他の「modality」の波形、複数回行なう再現性など用いて、相対的に、また、統計をもとに総合的に判断される。そこに調整、経験が必要となる（機械学習を適用することも可能である）。これは、幅、範囲を広げると他の機能の要素が入り始める一方、狭すぎると機能自体の要素がそぎ落とされてしまうので、そのレンジに関しては、調整が必要となるからである。例えば、複数回のデータがあると、レンジ、Hzと測定一致幅などが規定しやすい。
【0060】
［同調一致率について］
本明細書では、画像変化の傾向を、同調一致率として説明する。例えば、肺野を検出し、複数のブロック領域に分割し、各フレーム画像におけるブロック領域の「平均density（画素値ｘ）」を算出する。そして、「平均density（画素値ｘ）」の最小値から最大値の変化幅（０％～１００％）に対する各フレーム画像におけるブロック領域の平均画素値の割合（ｘ’）を算出する。一方、横隔膜の最小位置から最大位置の変化幅（０％～１００％）に対する各フレーム画像の横隔膜の変化（ｙ）の割合（ｙ’）との比の値（ｘ’／ｙ’）を用いて、比の値（ｘ’／ｙ’）が予め定められた一定の範囲内にあるブロック領域のみを抽出する。
本明細書では、画像変化の傾向を、同調一致率として説明する。例えば、肺野を検出し、複数のブロック領域に分割し、各フレーム画像におけるブロック領域の「平均density（画素値ｘ）」を算出する。そして、「平均density（画素値ｘ）」の最小値から最大値の変化幅（０％～１００％）に対する各フレーム画像におけるブロック領域の平均画素値の割合（ｘ’）を算出する。一方、横隔膜の最小位置から最大位置の変化幅（０％～１００％）に対する各フレーム画像の横隔膜の変化（ｙ）の割合（ｙ’）との比の値（ｘ’／ｙ’）を用いて、比の値（ｘ’／ｙ’）が予め定められた一定の範囲内にあるブロック領域のみを抽出する。
【0061】
ここで、ｙ’＝ｘ’若しくはｙ＝ａｘ（ａは横隔膜の振幅の数値や「density」の数値の係数）となる場合は、完全一致である。しかし、完全一致の場合のみが有意義な値であるわけではなく、ある一定の幅を持った値を抽出すべきである。そこで、本発明の一態様では、対数（log）を用いて、一定の幅を以下のように定める。すなわち、y=xの割合（%）で計算すると、同調の完全一致は「log y’/x’=0」である。さらに、同調一致率の範囲が狭いもの（数式的に狭い）範囲を抽出する場合は、例えば、０に近い範囲で「log y’/x’=-0.05~+0.05」と定め、同調一致率の範囲が広いもの（数式的に広い）範囲のものであれば例えば０に近い範囲で「log y’/x’=-0.5~＋0.5」と定める。すなわち、同調率の対数の値が、０を含む一定の範囲として定められる。この範囲が狭いものであればあるほど、また、その範囲内で一致する数値が高いほど、一致率が高いといえる。
【0062】
画素のピクセルごとにこの比の値を求めて個数をカウントすると、健康な人の場合は、完全一致の場合をピークとした正規分布が得られる。これに対し、疾患を有する人の場合は、この比の値の分布が崩れることとなる。なお、上記のように、対数を用いて幅を定める手法は、あくまでも一例であり、本発明はこれに限定されない。すなわち、本発明は、（ある大雑把な範囲の「density」の変化）≒（胸郭の変化）≒（横隔膜の動き）≒（肺機能検査）≒（胸腹呼吸センサの動き）≒（肺野の面積及び体積）として、“画像抽出”を行なうものであり、対数を用いる手法以外の手法も適用可能である。このような手法により、同調性画像を表示することが可能となる。
【0063】
血管の場合は、一連の心臓の収縮（ｙ）に呼応して生じる一連の「density」の変化（ｘ（肺門部における一波形）において、そのままの形でわずかな時間の遅れ（位相の変化）が存在するため、ｙ＝ａ’（ｘ－ｔ）と表される（すなわち、Y≒X）。完全一致の場合は、ｔ＝０であるため、ｙ＝ｘ、またはｙ＝ａ’ｘである。横隔膜の場合と同様に、同調一致率の範囲が狭いもの（数式的に狭い）範囲を抽出する場合は、例えば、０に近い範囲で「log y’/x’=-0.05~+0.05」と定め、同調一致率の範囲が広いもの（数式的に広い）範囲のであれば例えば０に近い範囲で「log y’/x’=-0.5~＋0.5」と定める。この範囲が狭いものであればあるほど、また、その範囲内で一致する数値が高いほど、一致率が高いといえる。
【0064】
その他の血管の場合は、上記の「心臓に呼応する部分」が除外され、肺門からプロットした中枢側の「density」を用いる。末梢の血管の場合も同様に取り扱うことができる。
【0065】
さらに、循環器についても本発明を適用することができ、例えば、心臓の「density」の変化が、肺門部～末梢肺野への血流の「density」の変化に直接関連し、一連の心臓の「density」の変化や肺門部の「density」の変化は、一種の変換を受けてそのまま伝播される。それは、心臓の「density」の変化と肺門部の「density」の変化の関係より若干の位相の差を得て生じると考えられる。また、肺門部などの「density」の変化が、そのまま肺野の血流への「density」の変化に関連するので、そのままの率で反映したもの（Y≒Xの一致率の関係）で同調性を表現することも可能である。また、頸部血管系や、胸部、腹部、骨盤、四肢などの大血管系においても、同様に、近傍の中枢の心臓血管でプロットした「density」の変化が、直接関連し、またはわずかな位相を伴って関連していると考えられる。そして、その「density」が、背景に応じて変動し、伝搬するときには「density」の変化の様が伝わるとして、同調一致率として考察することが可能となる。
【0066】
ここで、１枚の画像の変化量と１枚の画像の変化率のそれぞれにおいて、「吸気量合計≒呼気量合計」とすることができる。そこで、周囲空気の透過性との差から相対的に数値を出す場合、肺野の「density」から変化量を１としたときの相対的な値（Standard Differential Signal Density/Intensity）として表示しようとすると、（１）画像１枚ごとの差の画像で、１枚ごとに１とした時の画像（通常想定）、（２）１枚ごとの差の画像で「density（変化量や変化率）」を足した吸気全体もしくは呼気全体、もしくは吸気呼気の絶対値を１としたときの割合、さらに、（３）複数回の撮影におけるそれぞれの呼吸時（10%の時を数回select）における「density」総量を１としてその割合、について、それぞれ、変化量、変化率の描出を行なうことができる。
【0067】
また、ＭＲなどの３Ｄの場合であるが、吸気全体の「intensity（ＭＲの場合）」や「density」（ＣＴの場合）を合計した値（その際はそれを１としたとき）、その「intensity」や「density」の差は吸気（安静時や努力呼吸でも）の「peak flow volume deta」に換算でき、その値をその「intensity」や「density」の割合を出すことによって、少なくともＭＲＩやＣＴなどでの「３D×time」の計算する場合に、各肺野の部分における実測呼吸量、呼吸率を換算することができる。同様に、１回心拍出量を入力することで、肺野の「flow」における「capillary phase」における分布が、肺血流末梢量の分布、容量に換算する推定値を提示することも可能である。
【0068】
すなわち、（画像１枚あたりの吸気変化量）×（吸気の枚数全部）≒（画像１枚当たりの呼気変化量）×（呼気の枚数全部）≒（その時の吸気呼吸：自然呼吸もしくは努力呼吸のvolume）≒（その時の呼気呼吸：自然呼吸もしくは努力呼吸のvolume）≒（その時の自然呼吸もしくは努力呼吸の「volume」の吸気もしくは呼気の変化量）が成り立つ。１０％や２０％の１枚の変化量だけを取り出す場合は、（すべての枚数）×（その時間の変化量）を計算することによって、推定値を算出することが可能である。
【0069】
その抽出変化量を可視化して、画像に描出する。これが、以下に説明する呼吸機能解析、血管解析である。そして、胸郭や横隔膜の変化率を可視化する。その際に、再度結果に対するアーティファクトを除外し、新しいデータ抽出波形や最初のベースとなるデータ波形、その他のモダリティなどの波形、周囲、複数回の波形から抽出して、機能の抽出を行なう場合もある。アーティファクトを除外する手法は、後述する。
【0070】
また、上記抽出したもの以外から抽出した変化成分を除外したものでも特徴量を把握する場合がある。例えば、腹部腸管の動きを把握する際、腹部から呼吸の影響と血管の影響を除外して、腹部腸管の動きの抽出を図る。
【0071】
以下、本発明の実施形態について図面を参照して説明する。図１Ａは、本実施形態に係る診断支援システムの概略構成を示す図である。この診断支援システムは、コンピュータに診断支援プログラムを実行させることにより特定の機能を発揮する。基本モジュール１は、呼吸機能解析部３、肺血流解析部５、その他の血流解析部７、フーリエ解析部９、波形解析部１０および視覚化・数値化部１１から構成されている。基本モジュール１は、入力インタフェース１３を介してデータベース１５から画像データを取得する。データベース１５には、例えば、ＤＩＣＯＭ（Digital Imaging and COmmunication in Medicine）による画像が格納されている。基本モジュール１から出力された画像信号は、出力インタフェース１７を介してディスプレイ１９に表示される。次に、本実施形態に係る基本モジュールの機能について説明する。
【0072】
［呼吸要素の周期解析］
本実施形態では、以下の指標に基づいて呼吸要素の周期を解析する。「呼吸要素」とは、上述したように、呼気または吸気の全部または一部を含む概念である。すなわち、肺野内のある一定領域における「density」／「intensity」、横隔膜の動き、胸郭の動きの少なくとも一つを用いて呼吸要素の少なくとも一つの周波数を解析する。この「呼吸要素の少なくとも一つの周波数」において、呼吸要素が示す周波数スペクトルは一つ以上であり、一定の帯域幅を有する場合を含む概念である。肺野をブロックの集合体と考え、各ブロックから複数の周波数が抽出されることから、本実施形態では、これらを周波数群として処理する。なお、前述したように、ベースデータは、「波の形態自体」、および、「波の間隔（周波数：Hz）」の両方の概念を有するため、波の形態として処理することも可能である。また、Ｘ線（その他ＣＴ、ＭＲＩなどの複数種類のモダリティ）の透過性が高い部位で測定されるある一定のｖｏｌｕｍｅ 「density」／「intensity」で構成される範囲、スパイログラムなどの他の測定方法から得られるデータや外部入力情報を用いても良い。
本実施形態では、以下の指標に基づいて呼吸要素の周期を解析する。「呼吸要素」とは、上述したように、呼気または吸気の全部または一部を含む概念である。すなわち、肺野内のある一定領域における「density」／「intensity」、横隔膜の動き、胸郭の動きの少なくとも一つを用いて呼吸要素の少なくとも一つの周波数を解析する。この「呼吸要素の少なくとも一つの周波数」において、呼吸要素が示す周波数スペクトルは一つ以上であり、一定の帯域幅を有する場合を含む概念である。肺野をブロックの集合体と考え、各ブロックから複数の周波数が抽出されることから、本実施形態では、これらを周波数群として処理する。なお、前述したように、ベースデータは、「波の形態自体」、および、「波の間隔（周波数：Hz）」の両方の概念を有するため、波の形態として処理することも可能である。また、Ｘ線（その他ＣＴ、ＭＲＩなどの複数種類のモダリティ）の透過性が高い部位で測定されるある一定のｖｏｌｕｍｅ 「density」／「intensity」で構成される範囲、スパイログラムなどの他の測定方法から得られるデータや外部入力情報を用いても良い。
【0073】
なお、一呼吸毎の解析結果を比較し、複数のデータから傾向を解析して、データの確度を高めることもできる。
【0074】
また、呼吸要素の少なくとも一つの周波数の位相を変化させ、または呼吸要素の波形を円滑化させることで、呼吸要素を補正することも可能である。この場合、（胸郭、その他横隔膜の動き）≒（胸郭の動き）≒（density）≒（精密肺機能）≒（胸郭センサ）などの動きを用いてその波に位相を合わせる。また、肺野平均の「density」を追跡し、最後の変化は波の形態として波の二乗法などの近似を行ない、波の同定を行なう。ここで、胸部の「density」などの場合、最も変化する値が肺の「density」であるため、画面全体の「density」の評価により、肺の「density」の変化を評価する場合もある。波のプロットをする場合、実際動いている場合と、計測値で位相のずれが生じる場合がある。その場合は位相差を最大、最小値の位置、波の形態全体などで位相を補正する場合がある。
【0075】
［波形解析］
呼吸要素の波形から、波形の周波数の構成要素を算出することができる。これにより、上述した「波形同調性画像」を取得する。具体的には、解析範囲内のいずれかの部位の波形を特定し、前記特定した波形の周波数の構成要素を抽出し、前記波形の周波数の構成要素に対応する画像を出力する。
呼吸要素の波形から、波形の周波数の構成要素を算出することができる。これにより、上述した「波形同調性画像」を取得する。具体的には、解析範囲内のいずれかの部位の波形を特定し、前記特定した波形の周波数の構成要素を抽出し、前記波形の周波数の構成要素に対応する画像を出力する。
【0076】
［心血管拍解析および血管拍解析］
本実施形態では、以下の指標に基づいて心血管拍解析および血管拍を解析する。すなわち、心電図や脈拍計等の他のモダリティの計測結果、または肺輪郭から心臓・肺門位置・主要血管を特定し、各部位の「density」／「intensity」の変化を用いて血管拍を解析する。また、マニュアルで画像上にプロットし、対象部位の「density」／「intensity」の変化を解析しても良い。そして、心拍または血管拍から得られる心血管拍要素の少なくとも一つの周波数（波形）を特定する。なお、一拍毎の解析結果を比較し、複数のデータから傾向を解析して、データの確度を高めることが望ましい。また、各部位の「density」／「intensity」の抽出は、複数回実施したり、一定の範囲に対して行なうことで精度を高めることが可能となる。また、心血管拍周波数もしくは周波数帯を入力する方法もある。
本実施形態では、以下の指標に基づいて心血管拍解析および血管拍を解析する。すなわち、心電図や脈拍計等の他のモダリティの計測結果、または肺輪郭から心臓・肺門位置・主要血管を特定し、各部位の「density」／「intensity」の変化を用いて血管拍を解析する。また、マニュアルで画像上にプロットし、対象部位の「density」／「intensity」の変化を解析しても良い。そして、心拍または血管拍から得られる心血管拍要素の少なくとも一つの周波数（波形）を特定する。なお、一拍毎の解析結果を比較し、複数のデータから傾向を解析して、データの確度を高めることが望ましい。また、各部位の「density」／「intensity」の抽出は、複数回実施したり、一定の範囲に対して行なうことで精度を高めることが可能となる。また、心血管拍周波数もしくは周波数帯を入力する方法もある。
【0077】
［肺野同定］
データベース（DICOM）から画像を抽出し、上記の呼吸要素の周期解析結果を用いて、肺輪郭を自動検出する。この肺輪郭の自動検出については、従来から知られている技術を用いることができる。例えば、特開昭６３－２４０８３２号公報、または特開平２－２５０１８０号公報に開示されている技術を用いることが可能である。次に、肺野を複数のブロック領域に分けて、各ブロック領域の変化を計算する。ここで、撮影速度に応じてブロック領域の大きさを定めても良い。撮影速度が遅い場合は、あるフレーム画像の次のフレーム画像で対応する部位が特定しにくくなるため、ブロック領域を大きくする。一方、撮影速度が速い場合は、単位時間当たりのフレーム画像数が多いため、ブロック領域が小さくても追従することが可能となる。また、呼吸要素の周期のうちどのタイミングを選ぶかに応じて、ブロック領域の大きさを計算しても良い。ここで、肺野領域のずれを補正することが必要になる場合がある。その際には、胸郭の動き、横隔膜の動き、肺野全体の血管の位置関係を同定し、また、肺輪郭の相対位置を把握し、その動きに基づいて相対的に評価する。なお、ブロック領域が小さすぎると、画像のちらつきが発生する場合がある。これを防止するため、ブロック領域は一定の大きさを有する必要がある。
データベース（DICOM）から画像を抽出し、上記の呼吸要素の周期解析結果を用いて、肺輪郭を自動検出する。この肺輪郭の自動検出については、従来から知られている技術を用いることができる。例えば、特開昭６３－２４０８３２号公報、または特開平２－２５０１８０号公報に開示されている技術を用いることが可能である。次に、肺野を複数のブロック領域に分けて、各ブロック領域の変化を計算する。ここで、撮影速度に応じてブロック領域の大きさを定めても良い。撮影速度が遅い場合は、あるフレーム画像の次のフレーム画像で対応する部位が特定しにくくなるため、ブロック領域を大きくする。一方、撮影速度が速い場合は、単位時間当たりのフレーム画像数が多いため、ブロック領域が小さくても追従することが可能となる。また、呼吸要素の周期のうちどのタイミングを選ぶかに応じて、ブロック領域の大きさを計算しても良い。ここで、肺野領域のずれを補正することが必要になる場合がある。その際には、胸郭の動き、横隔膜の動き、肺野全体の血管の位置関係を同定し、また、肺輪郭の相対位置を把握し、その動きに基づいて相対的に評価する。なお、ブロック領域が小さすぎると、画像のちらつきが発生する場合がある。これを防止するため、ブロック領域は一定の大きさを有する必要がある。
【0078】
上記自動検出した肺野領域に少なくとも一つのベジエ曲線を用いて、肺野を点および制御点の座標として表すことができる。そして、少なくとも一つのベジエ曲線を用いて囲まれた閉じた曲線、いわゆる「単純閉曲線」を、複数用いることによって、肺野を表すことも可能である。同様に、一つまたは複数の単純閉曲線を用いて解析対象を表すことも可能である。
【0079】
各フレームの肺野は、特定のフレームにおいて検出した肺野上の少なくとも一つ以上のベジエ曲線（Bezier curve）を用いて、他のフレームにおける肺野を検出することもできる。例えば、最大と最小の２つの肺野を検出し、その２つの肺野を用いて、その他のフレームの肺野を算出する方法が挙げられる。ここでは、その他のフレームに「変化率」という変数を定義する。「変化率」は、肺野の大きさ、すなわち呼吸の状態を表現する値であり、横隔膜の位置や画像全体の「intensity」平均値などから算出できる。スパイログラフィーなどの外部装置の計測データを用いて算出したり、モデル化された変化率を用いたりすることも可能である。このように、「変化率」という変数を任意に定めることができるため、例えば、肺野が一定の割合（10%, 20%, 30%…）で変化していると仮定して算出することもできる。このように定義した変化率は、誤差を含んでいる場合があるため、誤差の自動・手動除去を行なった結果、または最小二乗法等を用いて近似を行なった結果などを使用して以後の処理を行なう場合もある。最大肺野と最小肺野まで、線形に変形すると仮定し、それぞれのフレームの変化率を用いて、線形変換などの技法を用いてそれぞれのフレームにおける肺野を算出する。
【0080】
また、連続するフレームの任意の範囲において、上述の処理を適用することが可能である。例えば、呼吸において、肺野は極大と極小への変化を繰り返すものであるが、実際の測定においては、極大時の形状が常に一定とは限らない。例えば、極大から極小、極小から極大の各範囲において、上述の処理を適用することによって、最大と最小の２つの肺野を定義して計算するよりも、精度高く肺野を算出することができると期待される。なお、ここでは、具体例として、極大と極小を用いて説明したが、本発明はこれに限定されるわけではなく、「任意の範囲」であるため、例えば、呼吸の途中、０％と３０％、３０％と１００％といった位置で行なうことも可能である。
【0081】
また、精度は低下するものの、１つの肺野から各フレームの肺野を算出することも可能である。例えば、胸郭の形状等から類推することで、肺野の変化ベクトルを規定する。具体的には、ベジエ曲線の制御点のそれぞれに変化ベクトルを規定する手法を採るが、本発明は、これに限定されるわけではない。そして、検出した１つの肺野と変化ベクトル、それぞれのフレームにおける変化率を使用して、それぞれのフレームにおける肺野を算出する。この算出結果に対して自動、または手動で補正を行なうことでさらに精度を高めることができる。また、３Ｄにおいても、上述の手法は有効である。すなわち、３Ｄの場合であっても、特定のフレームにおいて検出した肺野上の少なくとも一つ以上のベジエ曲面（Bezier surface）を用いて、他のフレームにおける肺野を検出する処理を実行することも可能である。これにより、フレーム間の肺野の画像を得ることが可能となる。
【0082】
図６Ｃは、呼吸要素の周期を示すグラフである。図６Ｃの画像中、白い垂直線が示されているが、これは、呼吸要素の周期中、現時点の位置を示す直線（指標）であり、図６Ｂに示す肺の動画の動きに応じて、呼吸要素の周期中の現在位置を示すように動く。呼吸要素の周期の現在位置を示すことによって、肺の動きの周期における現在位置を明確に把握することが可能となる。なお、本発明においては、呼吸要素の周期をグラフで表すだけではなく、血流の「density」、胸郭、横隔膜など、肺の動きと連動するものについては、すべてグラフ化することが可能である。
【0083】
また、被写体が「息を止めた場合」は、呼吸要素の周波数を特定できない場合がある。この場合は、被写体の心拍または血管拍から抽出される心血管拍要素の少なくとも一つの周波数を用いて、後述するフーリエ解析を行なう。この場合、心臓、横隔膜、または呼吸と連動する動的部位の動き方に応じて、後述するブロック領域の分割の仕方を変えるようにしても良い。
【0084】
［辺縁の検出とその評価］
本発明は肺の辺縁を検出し、その辺縁を評価することが可能である。例えば、前述の方法で肺野を算出した後、辺縁の位置および形状を、改めて精度高く検出することができる。算出された肺野内の任意の位置に点をプロットして、そこから四方八方に線を延伸し、各線において画素値の変化を評価する。例えば、図１４に示すように、肺を切断する線分Ｓに沿って画素値を算出すると、辺縁で画素が大きく変動することがわかるが、その変動の絶対値は異なる。例えば、左側の辺縁と右側の辺縁検出時の閾値を調整することで、辺縁検出の精度が高まる。また、領域ごとの画素値変動の特性を利用することもできる。図１４に示すように、Ｓ２領域とＳ３領域の縁の差分が小さかったとしても、画素値の変動の分散からＳ２領域とＳ３領域の縁を特定することができる。ここでは分散に着目したが、本発明はこれに限定されるわけではない。
本発明は肺の辺縁を検出し、その辺縁を評価することが可能である。例えば、前述の方法で肺野を算出した後、辺縁の位置および形状を、改めて精度高く検出することができる。算出された肺野内の任意の位置に点をプロットして、そこから四方八方に線を延伸し、各線において画素値の変化を評価する。例えば、図１４に示すように、肺を切断する線分Ｓに沿って画素値を算出すると、辺縁で画素が大きく変動することがわかるが、その変動の絶対値は異なる。例えば、左側の辺縁と右側の辺縁検出時の閾値を調整することで、辺縁検出の精度が高まる。また、領域ごとの画素値変動の特性を利用することもできる。図１４に示すように、Ｓ２領域とＳ３領域の縁の差分が小さかったとしても、画素値の変動の分散からＳ２領域とＳ３領域の縁を特定することができる。ここでは分散に着目したが、本発明はこれに限定されるわけではない。
【0085】
さらに、同様の考え方によって、肺以外の臓器、血管、腫瘍などの解析範囲の辺縁を検出することが可能になる。例えば、血管中に造影剤が存在する場合、血管の内部は明確に可視化できるが、血管の外側や厚さを明確に算出することは容易でない。本実施形態では辺縁を正確に検出できるため、解析範囲内にある血管の形状、特徴を算出することができる。これにより、従来は把握が容易でなかった血管の厚さや外周を定量的にも把握し、診断に用いることが可能となる。
【0086】
［ブロック領域の作成］
肺野を複数のブロック領域に分ける手法について説明する。図１Ｂは、肺野を肺門から放射状に分割する手法を示す図である。肺は、肺尖側よりも横隔膜側の方が大きく動くため、横隔膜側に近いほど粗く分割した点をプロットするようにしても良い。なお、図１Ｂにおいて、縦方向の線（点線）を追加的に描画し、複数の矩形（正方形）のブロック領域に分けてもよい。これにより、肺の動作をより正確に表すことが可能となる。なお、「肺の縦方向に点をプロットし、肺を横断的に分ける手法」、「肺の横方向に点をプロットし、肺を縦断的に分ける手法」、「肺尖部における接線と横隔膜における接線を引き、その接線が交わる点を中心点として定め、その点を含む直線（例えば、鉛直線）からある一定の角度で引いた線分で肺を分割する手法」、「肺を肺尖（または肺門）から横隔膜端部を結ぶ直線と直交する複数の平面で切断する手法」などの手法で、肺野を分割することも可能である。なお、これらの手法は、三次元立体画像にも適用可能である。３Ｄの場合は、複数の曲面または平面で囲まれた空間として、各臓器を捉える。臓器をさらに細かく分けることもできる。例えば、右肺の３Ｄモデルを考えた場合、上葉、中葉、下葉に分けて取り扱うこともできる。
肺野を複数のブロック領域に分ける手法について説明する。図１Ｂは、肺野を肺門から放射状に分割する手法を示す図である。肺は、肺尖側よりも横隔膜側の方が大きく動くため、横隔膜側に近いほど粗く分割した点をプロットするようにしても良い。なお、図１Ｂにおいて、縦方向の線（点線）を追加的に描画し、複数の矩形（正方形）のブロック領域に分けてもよい。これにより、肺の動作をより正確に表すことが可能となる。なお、「肺の縦方向に点をプロットし、肺を横断的に分ける手法」、「肺の横方向に点をプロットし、肺を縦断的に分ける手法」、「肺尖部における接線と横隔膜における接線を引き、その接線が交わる点を中心点として定め、その点を含む直線（例えば、鉛直線）からある一定の角度で引いた線分で肺を分割する手法」、「肺を肺尖（または肺門）から横隔膜端部を結ぶ直線と直交する複数の平面で切断する手法」などの手法で、肺野を分割することも可能である。なお、これらの手法は、三次元立体画像にも適用可能である。３Ｄの場合は、複数の曲面または平面で囲まれた空間として、各臓器を捉える。臓器をさらに細かく分けることもできる。例えば、右肺の３Ｄモデルを考えた場合、上葉、中葉、下葉に分けて取り扱うこともできる。
【0087】
肺野領域は、胸郭の動き、横隔膜の動き、肺野全体の血管の位置関係を同定し、肺輪郭の相対的位置を把握し、その動きに基づいて相対的に評価すべきものである。このため、本願発明では、肺輪郭を自動検出した後、肺輪郭によって特定される領域を複数のブロック領域に分割し、各ブロック領域に含まれる画像の変化の値（画素値）を平均化する。例えば、図１０に示すように、ベジエ曲線上で対向する肺の辺縁上に点をプロットし、それらを接続した上で、その中間の点を通る曲線を用いることも可能である。その結果、図１Ｃに示されているように、肺の形態が時間の経過によって変化しても、注目する領域の経時的変化を追跡することが可能となる。一方、図１Ｄは、解析対象となる臓器（この場合は肺）の形態を考慮せずに、ブロック領域に分割した場合の経時的変化を示す図である。肺野領域とは、上述した通り、胸郭の動き、横隔膜の動き、肺野全体の血管の位置関係を同定し、肺輪郭の相対的位置を把握し、その動きに基づいて相対的に評価すべきものであるが、図１Ｄに示すように、肺野領域を特定することなくブロック領域に分割すると、肺の経時的変化により、注目領域が肺野領域から外れ、意味のない画像となってしまう。特に、横隔膜の動きは肺野を縮める動作が強いため、横隔膜や全体的な数値だけを補正するのではなく、胸郭成分やその他の複数の要素を取り込んで、肺野領域を補正することが好ましい。また、呼吸要素周波数もしくは周波数帯を入力する方法もある。３Ｄについても、同様に領域分割の計算ができる。
【0088】
さらに、図１１に示すように、肺野Ａにおいて、ベジエ曲線を用い、検出した肺野内に内部制御点を選定し、肺野内の内部制御点を通る曲線または直線によって肺野を分割することも可能である。すなわち、肺野の枠だけではなく、肺野領域の内部にも制御点を設け、これらの制御点を用いて肺野領域（Ａ）を分割する。この場合、図１２に示すように、検出した肺野の外延およびその近傍における制御点の間隔を相対的に大きくし（１）、検出した肺野内における部位毎の膨張比率に応じて、内部制御点の間隔を相対的に小さく（２）するようにしても良い。また、肺野Ａにおいて、制御点の間隔を、人体に対して頭尾方向に進むに従って相対的に大きくしても良く、または、特定のベクトル方向に従って相対的に大きくしても良い。このベクトルの定め方は任意であるが、例えば、肺尖から肺野の反対側へ向かう方向に定めても良いし、図１Ｂに示すように、肺門から肺野の反対側へ向かう方向に定めても良い。また、肺の構造に応じた方向にベクトルを定めることも可能である。このように、肺野の分割の仕方を「不等分割」とすることによって、領域毎の特徴を考慮した画像の表示が可能となる。例えば、肺野の外周は動きが大きく、ずれが大きくなるため、ブロックを大きくする一方、肺野の内部は動きが小さく、ずれも小さくなるため、ブロックを小さくして細かくする。また、例えば、肺野の横隔膜側は動きが大きく、ずれが大きくなるため、ブロックを大きくする一方、肺野の頭部側は動きが小さく、ずれも小さくなるため、ブロックを小さくして細かくする。これにより、表示の精度を高めることが可能となる。この手法は、肺野に限定されるわけではなく、呼吸と連動する動的部位などにも適用することが可能である。このような手法は、肺を肺葉ごとに３次元的に分割する場合にも適用できる。また、横隔膜の下側の部位、例えば、心臓やその他の臓器をベジエ曲線で囲って表示する場合に用いることも可能である。この場合も、心臓やその他の臓器の構造に応じた方向にベクトルを定めて、領域を不等分割することも可能である。
【0089】
次に、アーティファクトを排除して画像データを補間する。すなわち、解析範囲内に骨などが含まれるとノイズとして表れてしまうため、ノイズカットフィルタを用いてノイズを除去することが望ましい。Ｘ線画像においては、通例では、空気を－１０００とし、骨を１０００としているため、透過性が高い部分は画素値が低く、黒く表示され、透過性の低い部分は、画素値が高く、白く表示される。例えば、画素値を２５６階調で表す場合、黒は０で白は２５５となる。肺野領域内で、血管や骨が存在しない位置の周辺は、Ｘ線が透過しやすいため、Ｘ線画像の画素値が低くなり、Ｘ線画像は黒くなる。一方、血管や骨が存在する位置は、Ｘ線が透過し難いため、Ｘ線画像の画素値が高くなり、Ｘ線画像は白くなる。その他ＣＴ、ＭＲＩにおいても同様のことが言える。ここで、上記の呼吸要素の周期解析の結果から、一呼吸あたりの波形に基づいて、同一位相の値を用いてデータを補間し、アーティファクトを排除することが可能となる。また、「座標が異なること」、「画素値が極端に変動すること」、「周波数やdensityが異常に高くなること」を検出した場合に、それらに対してカットオフを行ない、残りの得られた画像に対して、例えば、最小二乗法等を用いて連続的で滑らかな波の形を同定することによって、横隔膜のHz計算、肺野の調節に使用できるようにしても良い。また、画像を重ね合わせる場合、（１）前後で片方の画像を取得した取得比較画像をその座標のまま重ねる場合と、（２）前後片方の画像をbaseに取得した後、画像を相対的に拡張してその相対的位置情報をbaseに重ねる方法とがある。以上のような手法によって、肺野の形態を修正したり、ブロック領域の画像の変化を修正したりすることが可能となる。その際、再度、結果に対するアーティファクト（artifact）を除外し、新しいデータ抽出波形や最初のベースとなるデータ波形、その他のモダリティなどの波形、周囲、複数回の波形から抽出して、機能の抽出を行なう。その際、回数は一回でも複数回でも良い。
【0090】
ここで、時間軸における「reconstruction」について説明する。例えば、１５ｆ/ｓの吸気時間が２秒の場合、３０＋１枚の画像が得られる。その場合、１０％ずつの「reconstruction」は、単純に3枚ずつ重ね合わせれば実施できる。その際、例えば、０．１秒が１０％で、その画像が０．０７秒と０．１２秒の写真しか取得していない場合は、０．１秒の「reconstruction」が必要となる。その場合、１０％前後の画像の中間の値（両者の平均）値を与えて「reconstruction」を行なう。また、時間軸で捉え、その時間の割合で係数を変えても良い。例えば、時間軸の差があって、０．１秒の撮影の値がなく、０．０７秒と０．１２秒の撮影時間があるときは、「（その0.07秒の値）×2/5＋（0.12秒の値）×3/5」と計算し直して、「reconstruction」を行なうことができる。さらに、呼吸の平均や横隔膜の係数の変化量からその秒における変化位置関係を認識し、その値を係数にして数字の割合を求める。なお、「Maximum Differential Intensity Projection」の０～１００％を含み、１０％から２０％の「reconstruction」や、１０％から４０％の「reconstruction」など、厚みを持たせて計算することが望ましい。このように、撮影していない部分についても、１呼吸割合での「reconstruction」を行なうことが可能である。なお、本発明は、呼吸のみならず、血流、胸郭の動き、横隔膜、その他これらと連動する一連の動きに対しても同様に「reconstruction」を行なうことが可能である。また、ブロックごと、またはピクセルごとに「reconstruction」をすることも可能である。なお、「Maximum Differential Intensity Projection」の０～１００％を含み、１０％から２０％の「reconstruction」や、１０％から４０％の「reconstruction」など、厚みを持たせて計算することが望ましい。
【0091】
また、上述した手法で肺野を検出し、この検出した肺野を正規化しても良い。すなわち、検出した肺野を空間的に正規化し、またはリコンストラクション（reconstruction）を利用して時間的に正規化する。人体の相違によって肺野の大きさや形状が異なるが、これを正規化することで一定の領域内に表示することができる。
【0092】
［横隔膜および胸郭］
上記のように肺野を同定すると、横隔膜の動きや胸郭についても把握することが可能となる。すなわち、認識した横隔膜のＸｐ上（２Ｄ画像）の横隔膜による曲線や胸郭の曲線を細かな座標の集まりとして計算し、その平均や曲線の局部における下方への変化率や変化量、また横隔膜を曲線として「curve fitting」してその変形率を数値化することにより画像からの機能評価の位置づけを行なうことができる。また、横隔膜面以外の胸部で描いた辺縁の曲線についても、同様に細かな座標の集まりとして計算し、その平均や曲線の変化率を数値化することにより画像からの機能評価を行なうことができる。上記の２つの変化率、変化を、相対的・相互連動として評価し、変化率が異なる（同じように連動して動かない部位など）を数値化、画像化することによりｍｏｖｅｍｅｎｔの機能評価を行なう。
上記のように肺野を同定すると、横隔膜の動きや胸郭についても把握することが可能となる。すなわち、認識した横隔膜のＸｐ上（２Ｄ画像）の横隔膜による曲線や胸郭の曲線を細かな座標の集まりとして計算し、その平均や曲線の局部における下方への変化率や変化量、また横隔膜を曲線として「curve fitting」してその変形率を数値化することにより画像からの機能評価の位置づけを行なうことができる。また、横隔膜面以外の胸部で描いた辺縁の曲線についても、同様に細かな座標の集まりとして計算し、その平均や曲線の変化率を数値化することにより画像からの機能評価を行なうことができる。上記の２つの変化率、変化を、相対的・相互連動として評価し、変化率が異なる（同じように連動して動かない部位など）を数値化、画像化することによりｍｏｖｅｍｅｎｔの機能評価を行なう。
【0093】
ここで、「横隔膜および胸郭評価方法」について説明する。まず、横隔膜について、身体の軸（いわゆる正中線）に直交する左右水平線を軸にその動きを表示する。次に、横隔膜のラインを基線に平坦化する。すなわち、横隔膜のラインを水平の直線に合わせる。そして、横隔膜の動きを評価する。さらに、横隔膜のラインを伸ばして平坦化し、曲線の直交する動きを評価する。次に、胸郭外側では、肺尖から横隔膜胸郭角を結んだ線を基線として（軸として）動きを評価する。胸郭のラインを基線として平坦化し、すなわち、胸郭のラインを「肺尖－肋横隔膜角」との直線に合わせて、動きを評価する。胸郭のラインを基線に伸ばして平坦化して曲線の直交する動きを評価する。そして、上記胸郭、横隔膜ラインの曲率や曲率半径を評価する。そして、上記の変化を「変化量」として算出し、この変化量を微分して“変化率”として評価する。
【0094】
図６Ｂおよび図６Ｃは、ディスプレイに表示される画像の一例を示す図である。図６Ｂでは、左肺の動きが動画として表示される。図６Ｂの画像中、白い水平線が示されているが、これは、横隔膜の位置を示す直線（指標）であり、動画が再生されると、横隔膜の動きに追従して上下に動くこととなる。このように、横隔膜を検出し、検出した横隔膜の位置を示す指標、すなわち、横隔膜の位置を示す白い水平線を示すことによって、医師による画像診断を行なうことができるようになる。また、横隔膜の一部のみならず、肺野－横隔膜のラインの認識を使用し、すべての点を認識して、左右、外内側などの横隔膜の一領域、また、横隔膜全体の診断を行なうことが可能となる。同様に、横隔膜のみならず、呼吸と連動する動的部位、例えば、胸郭などの動きも同様に、接線位置などの直線や肺野認識による胸郭の直線によって、胸郭の動きを判断することが可能となる。このように、辺縁は動くものという仮定で、連続画像で差分を取ることによって、辺縁を検出することが可能となる。例えば、腫瘍は堅く、その周りは柔らかいことが多い。このため、腫瘍はあまり動かず、その周りは活発に動くこととなるため、差分を取ることによって、腫瘍の辺縁を検出することができる。
【0095】
また、ＭＲＩやＣＴなどの３Ｄ画像においても、横隔膜の面を一つの座標や立体的な曲面として捉え、その座標や曲面を細かな座標の集まり（横隔膜の辺縁の輪郭、平面および座標の集合群）として計算し、その平均や曲面の局部における下方への変化率や変化量、また横隔膜を曲面として「curve fitting」してその変形率を数値化することにより画像からの機能評価の位置づけを行なうことができる。また、横隔膜面以外の胸部で描いた辺縁の曲面についても、同様に細かな座標の集まりとして計算し、その平均や曲面の変化率を数値化することにより、画像からの機能評価を行なうことができる。上記２つの変化率、変化を相対的、相互連動として評価し、変化率が異なる（同じように連動して動かない部位など）を数値化、画像化することによりｍｏｖｅｍｅｎｔの機能評価を行なう。
【0096】
［フーリエ解析］
上記のように解析した呼吸要素の周期および血管拍周期に基づいて、各ブロック領域の「density」／「intensity」の値や、また、その変化量について、フーリエ解析を実施する。図２Ａは、特定ブロックの「intensity」変化と、それをフーリエ解析した結果を示す図である。図２Ｂは、心拍に近い周波数成分を抜き出したフーリエ変換結果と、これを逆フーリエ変換して心拍に近い周波数成分の「intensity」変化を示す図である。例えば、特定ブロックの「intensity」変化をフーリエ変換（フーリエ解析）すると、図２Ａに示すような結果が得られる。そして、図２Ａに示した周波数成分から、心拍に近い周波数成分を抜き出すと、図２Ｂの紙面に対して右側に示すような結果が得られる。これを逆フーリエ変換することによって、図２Ｂの紙面に対して左側に示すように、心拍の変化に同調した「intensity」変化を得ることができる。
上記のように解析した呼吸要素の周期および血管拍周期に基づいて、各ブロック領域の「density」／「intensity」の値や、また、その変化量について、フーリエ解析を実施する。図２Ａは、特定ブロックの「intensity」変化と、それをフーリエ解析した結果を示す図である。図２Ｂは、心拍に近い周波数成分を抜き出したフーリエ変換結果と、これを逆フーリエ変換して心拍に近い周波数成分の「intensity」変化を示す図である。例えば、特定ブロックの「intensity」変化をフーリエ変換（フーリエ解析）すると、図２Ａに示すような結果が得られる。そして、図２Ａに示した周波数成分から、心拍に近い周波数成分を抜き出すと、図２Ｂの紙面に対して右側に示すような結果が得られる。これを逆フーリエ変換することによって、図２Ｂの紙面に対して左側に示すように、心拍の変化に同調した「intensity」変化を得ることができる。
【0097】
図９に示すように、特定のスペクトルに係数をかけて重みを付けることも可能である。例えば、波形同調性を実現するために、この手法を用いることが可能である。すなわち、フーリエ逆変換を行なう時の周波数の選び方として、複数の周波数を選択して、かつその比率を掛けた後、フーリエ逆変換を行なう。例えば、抽出する帯域中、最も周波数の高いスペクトルを強調表示したい場合に、そのスペクトル強度を２倍することが可能となる。この場合、周波数の連続性は無くても構わない。とびとびに存在するスペクトルを選択することが可能である。
【0098】
また、心臓の「density」の位置を、左肺（内臓逆位などの場合右芯の場合もある）の形態（肺野抽出の形態から左肺のくぼんだ部位の領域）および椎体、横隔膜の位置から類推することが可能となる。この場合、心臓のROIをとって「density」の抽出を行なう。この抽出をする際には、呼吸、血流の相対的なスペクトル値でおおまかな領域を用いて類推する。また、予め、心血管拍で生じるＨｚ帯（心拍40～150Hz、≒0.67Hz～2.5Hz）などを用いて「filtering」を行なうことによって、呼吸やその他の「artifact」による周波数を除去する場合もある。また、心臓の位置も呼吸状況に応じて変わるので、胸郭の位置が変わるにつれて胸郭の形態値から相対的に心臓の位置を変更し、より正確な心血管拍の抽出や肺門、大血管などの抽出を行なうことがある。さらに、横隔膜の動きと同様、規則的に動く心臓の輪郭に基づいて、周波数を計算する方法がある。
【0099】
ここで、周波数成分からなるスペクトルに対して逆フーリエ変換を行なう際に、呼吸や血流の「density」から特定される周波数要素（呼吸周波数、心血管拍周波数）と、スペクトルの帯域（ＢＰＦ：band pass filterを用いても良い）とを両方を加味し、または、そのどちらかの要素に基づいて逆フーリエ変換を行なうようにしても良い。また、上記フーリエ変換後に得られるスペクトルから、臓器特有の周期的な変化のスペクトル構成比に基づいて、逆フーリエ変換を行なう際の少なくとも一つの周波数を選択してもよい。さらに、フーリエ変換後に得られる複数の周波数の構成割合によって、特定の臓器や解析対象となる領域の波形を特定することも可能である（波形同調性画像の作成）。
【0100】
なお、フーリエ変換を実行する際には、短時間で計算ができるように、ＡＲ法（Autoregressive Moving average model）を用いることが可能である。ＡＲ法では、自己回帰移動平均モデルにおいて、ユールウォーカー方程式（Yule-walker equiation）やカルマンフィルタを用いる方法があり、そこで導きだされるユールウォーカー推定値（Yule-walker estimates）、PARCOR法、最小二乗法を用いて、計算を補足することができる。これにより、より早く、リアルタイムに近い画像を取得したり、計算の補助やアーティファクト（artifact）の補正を行なったりすることが可能となる。このようなフーリエ解析により、各ブロック領域における画像の性質を抜き出して表示することが可能となる。
【0101】
また、このフーリエ解析の際に、「デジタルフィルタ」を用いる手法を採ることも可能である。すなわち、元の波形にフーリエ変換を行なって、各スペクトルのパラメータを取得し、元の波に演算処理を施す「デジタルフィルタ」を用いる。この場合は、逆フーリエ変換を行わず、デジタルフィルタを用いる。
【0102】
ここで、各フレーム画像における各ブロック領域の画像の変化をフーリエ変換し、フーリエ変換後に得られるスペクトルのうち、呼吸要素の周期に対応するスペクトルを含む一定の帯域内のスペクトルを抽出することができる。図２Ｃは、フーリエ変換後に得られたスペクトルのうち、ある一定の帯域を抽出する例を示す図である。合成波のスペクトルの周波数ｆは、合成元となる各周波数ｆ１（呼吸成分）、ｆ２（血流成分）との間に、「１／ｆ＝１／ｆ１＋１／ｆ２」という関係が成り立っており、スペクトルを抽出する際に、以下の方法を採ることが可能である。
【0103】
（１）血流のスペクトル比率が高い部分を抽出する。
（２）呼吸／血流に対応するスペクトルのピークとその近辺の複数の合成波のピークの中間で区切り、スペクトルを抽出する。
（３）呼吸／血流に対応するスペクトルのピークとその近辺の複数の合成波のスペクトルの谷の部分で区切り、スペクトルを抽出する。
（２）呼吸／血流に対応するスペクトルのピークとその近辺の複数の合成波のピークの中間で区切り、スペクトルを抽出する。
（３）呼吸／血流に対応するスペクトルのピークとその近辺の複数の合成波のスペクトルの谷の部分で区切り、スペクトルを抽出する。
【0104】
上述したように、本発明では、固定的なＢＰＦを用いているわけではなく、呼吸要素の周期に対応するスペクトルを含む一定の帯域内のスペクトルを抽出する。さらに、本願発明では、フーリエ変換後に得られるスペクトルのうち、フレーム画像から得られる呼吸要素以外の周波数（例えば、また、各部位の「density」／「intensity」、心拍または血管拍から得られる心拍要素）、またはオペレータによって外部から入力された周波数に対応するスペクトル（例えば、スペクトルモデル）を含む一定の帯域内のスペクトルを抽出することも可能である。
【0105】
ここで、合成波のスペクトルの成分は２つの成分（呼吸、血流）のみであれば５０％＋５０％となり、３つの成分の場合は、１／３ずつの配分となる。このため、呼吸成分のスペクトルが何％、血流成分のスペクトルが何％と、スペクトルの成分およびその高さからある程度、合成波のスペクトルを計算することができる。その割合（％）が高いところでスペクトルを抽出することが可能である。すなわち、血流成分／呼吸成分と合成波成分との割合を計算し、血流成分／呼吸成分の高いスペクトル値を計算して抽出する。なお、横隔膜の同定等において、呼吸や心臓血管の周波数を取得したdataから、Hz(周波数)が比較的一定になる部位、すなわち、Hzの変化が少ない領域に対応するスペクトルやその重ね合わせのみを抽出する場合もある。また、スペクトルの帯域を定める場合、横隔膜の同定等をする際に、Ｈｚの変化が生じたrangeおよびその周囲の領域でスペクトルの帯域を定める場合もある。波形の構成要素を加味することもある。
【0106】
なお、逆フーリエ変換する際のスペクトルについては、「単純モデル化した周波数および周波数帯から高い部位（一つもしくは複数）を用いて抽出する場合（シュミレーション主義）」と、「実際の周波数または周波数帯に基づいて、スペクトル値に応じて周波数の高い部位または周波数の低い部位を抽出する場合（現場主義）」とを選択することが可能である。また、心臓の周波数がＡ、肺の周波数がＢである場合、周波数帯域全体からＡを減算することによって、Ｂが得られる。また、フーリエ変換から取得されるスペクトルについては、周波数軸上の一か所のみならず複数の箇所を抽出することも可能である。
【0107】
以上により、呼吸要素の周期や血管拍周期に完全一致する場合のみならず、考慮した方が良いスペクトルも抽出することができ、画像診断に寄与することが可能となる。なお、「呼吸」や「心拍」は、特定の周波数帯域に含まれることが知られている。このため、呼吸の場合は、例えば、「0~0.5Hz（呼吸数0~30回/分）」、循環器の場合は、例えば、「0.6~2.5（心拍/脈拍数36~150回/分）Hz」というフィルタを用いて、予めこのフィルタで呼吸周波数や循環器の周波数を特定していくことも可能である。これにより、周波数同調性画像を表示することが可能となる。これは、心臓の「density」変化を取得する際に、呼吸（肺）の「density」変化を拾ってしまったり、肺の「density」変化を取得する際に、心臓の「density」変化を拾ってしまったりする場合があるからである。
【0108】
［視覚化・数値化］
上記のように解析した結果を、視覚化および数値化する。視覚化および数値化をする際に、本明細書では、「モデル化した肺」を定義する。肺を動画像で表示する際、位置関係が動いてしまうため、相対的判断が容易ではない。このため、位置関係のずれを、空間的に統一化・平均化する。例えば、肺の形状を扇形などの図形に当てはめ、形を整えた状態で表示する。そして、リコンストラクションの概念を用いて時間的に統一化する。例えば、「複数の呼吸のうち、20％の肺の状況」を抽出し、それを「一呼吸の20％の肺の状況」として定めることができる。このように、空間的、時間的に統一化した肺を「モデル化した肺」とする。これにより、異なる患者同士を比較したり、一人の患者の現在と過去とを比較したりする際に、相対的判断が容易となる。
上記のように解析した結果を、視覚化および数値化する。視覚化および数値化をする際に、本明細書では、「モデル化した肺」を定義する。肺を動画像で表示する際、位置関係が動いてしまうため、相対的判断が容易ではない。このため、位置関係のずれを、空間的に統一化・平均化する。例えば、肺の形状を扇形などの図形に当てはめ、形を整えた状態で表示する。そして、リコンストラクションの概念を用いて時間的に統一化する。例えば、「複数の呼吸のうち、20％の肺の状況」を抽出し、それを「一呼吸の20％の肺の状況」として定めることができる。このように、空間的、時間的に統一化した肺を「モデル化した肺」とする。これにより、異なる患者同士を比較したり、一人の患者の現在と過去とを比較したりする際に、相対的判断が容易となる。
【0109】
例えば、ｓｔａｎｄａｒｄ ｕｐｔａｋｅとして、計測された肺野全域の「density」／「intensity」から平均値を１として相対的／対数的に値を表示することがある。また、血流の方向だけを採用するため、特定方向への変化を切り出すことがある。これにより、意味のある方法のデータだけを取り出すことが可能となる。肺野同定結果を用いて、解析範囲の変化に追従して疑似カラー化を行なう。すなわち、フェーズに合わせた特定の形（最小、最大、平均、中央値）に沿って、各個人（被写体）の解析結果を相対的な領域に当てはめる。
【0110】
また、複数の解析結果を比較できる特定の形状・フェーズに変形させる。さらに、モデル化した肺を作成する際、上記呼吸要素の周期解析の結果を用いて、肺野内の相対的な位置関係を計算する。なお、モデル化した肺は、複数の患者の胸郭ライン、「density」、横隔膜などを総合的に平均化したラインを用いて作成する。モデル化した肺の作成の際、肺血流の場合は、肺門から肺端部にかけて放射状に距離を測ることができる。また、呼吸の場合は、胸郭や横隔膜の動きに応じて補正する必要がある。さらに、肺尖からの距離を考慮して複合的に計算しても良い。
【0111】
また、逆フーリエ変換後に、振幅値が相対的に大きいブロックのみを抽出して表示しても良い。すなわち、ブロック毎にフーリエ解析を行なう場合、逆フーリエ変換の後、波の振幅が大きいブロックと、波の振幅が小さいブロックが存在する。そこで、振幅が相対的に大きいブロックのみを抽出し、視覚化することも有効である。また、逆フーリエ変換後、各数値の実部と虚部とをそれぞれ使い分けることができる。例えば、実部のみから画像を再構成したり、虚部のみから画像を再構成したり、実部および虚部の絶対値から画像を再構成することが可能である。
【0112】
モデル化した肺に対して、フーリエ解析をしても良い。すなわち、呼吸数回の画像を合わせたり、フーリエ解析や相対位置把握をしたりする際にもモデル化した肺を用いることが可能である。モデル化した肺を用いることによって、取得した複数のフレームをモデル化した肺に当てはめたり、血管の場合、心拍（例えば、肺門部から得られる心拍など）に応じて計算したモデル化した肺に当てはめたりすることによって、フーリエ解析を行なう場合の相対位置を一定とすることが可能となる。ベースとなる呼吸の状態を取得する際に、モデル化した肺を使うことによって、安定的に計算結果を得ることも可能となる。また、肺をモデル化することで、空間の相違を固定化することができ、肺の動きを見やすくさせることが可能となる。
【0113】
画像化において、相対評価の標識方法は、以下の通りである。すなわち、画像を相対的に白黒、カラーマッピングで標識する。差分によって得られた「density」／「intensity」の数％前後の値をカットし、その上下残りを相対的に表示することがある。もしくは、得られた差分の前後数％前後の値はとびぬけた値となる場合があるため、これを「ａｒｔｉｆａｃｔ」として除去し、残りの部分を相対的に表示することがある。０～２５５諧調などの方法の他、０～１００％の値として表示することもある。
【0114】
なお、ピクセルをある程度曖昧に表示し、ぼやけた状態にして全体を表示することも可能である。特に、肺血管の場合、高い信号値の間に低い信号の値が混在するが、高い信号値のみを大雑把に把握することができれば、全体として曖昧であっても構わない。例えば、血流の場合は、閾値以上の信号を抜き出し、呼吸の場合は閾値以上の信号を抜き出さなかったりしてもよい。具体的には、次の表の数字を１ピクセルとして真ん中の数値を取得する場合、真ん中の数値が占める割合を取得し、１ピクセル内で平均化すると、隣接する画素との間で滑らかに表現することができる。この手法は、ブロックごとの平均インテンシティ算出の際にも用いることが可能である。
【表1】
【0115】
この手法は、肺野のみならず、任意の解析範囲のデンシティ（density）を検出し、デンシティが相対的に大きく変化する箇所を除去する際にも適用することが可能である。また、予め設定した閾値を、大きく超える点をカットオフする。また、肋骨の形態認識、例えば、突然出現する高／低信号ラインを認識し、除去する。また、位相から同様に、突然出現する信号、例えば、reconstructionのphaseが15%～20％前後でartifactが認められる患者の特徴など、通常の波の変化と異なる突然の信号を除去する場合がある。なお、最初にベースデータを取る際、（横隔膜）≒（胸郭）≒（胸郭の動き）≒（スパイロメータ）≒（肺野）、fieldの（density）≒（volumetry）などの計算において、位相が異なる場合があり、その位相を実際認識できる形態（XPの輪郭）に当てはめる場合がある。
【0116】
モデル化した肺が作成できると、上述したように、同調性、一致率、不一致率を数値化して提示することが可能となる（周波数同調性画像または波長同調性画像の表示）。これにより、正常な状態からの逸脱が表示できる。本実施形態によれば、フーリエ解析を実行することにより、新しい病気の可能性の発見、普通の自分との比較、手と足との比較や、反対側の手および足との比較が可能となる。さらに、足の動かし方、嚥下などでどこがおかしいのかを同調性の数値化で把握することが可能となる。また、病気の状態の人が一定時間経過後に変化しているかどうかを判断したり、変化している場合は、変化の前後を比較することが可能となる。また、肺野を、末梢からの距離を一定として放射状に見やすいような形態（円形～類円形）とすることで、内層～中層、外層などの評価をしやすくすることもでき、また、「血管の末梢優位」か「中層優位」かに応じて表現することも可能である。
【0117】
なお、視覚化の際に、フーリエ変換後の画像とフーリエ変換前の画像を切り替えて表示したり、両者を一つの画面に並べて表示したりすることも可能である。
【0118】
図２Ｄに示すように、モデル化した肺を１００としたときに、当該人体において、何パーセントの相違があるのかどうかを把握し、変化率を表示することが可能となる。なお、肺全体のみならず、肺の一部分でも相違を把握することが可能である。特に、上述したように、横隔膜の動きのみを特定すると共に、横隔膜以外の肺野の形状を固定して、横隔膜の動きを表示すると共に、同調一致率や変化率を表示することが可能である。さらに、肺野の全部を固定して、同調一致率や変化率を表示することも可能である。なお、「Variation分類」をすることによって、標準血流の特定も可能である。すなわち、呼吸要素の周期を特定し、血管の相対的な位置関係を計算し、被写体の血流動態を標準血流として特定することが可能となる。
【0119】
また、パターンマッチングの手法を用いて肺の検出をしても良い。図２Ｅ～図２Ｈは、肺野領域のパターン画像の例を示す図である。図２Ｅ～図２Ｈに示すように、肺の形を、パターン分類をしておき、これらのなかで近いものを抽出するようにしても良い。この手法により、対象の画像が、片肺を表すのか、両肺を表すのかを特定することができる。また、右肺であるのか、左肺であるのかを特定することもできる。パターン数は、限定しないが、４～５パターンを持っておくことが想定される。なお、このように、肺野の形態（形状）のみで右肺、左肺、両肺の認識をする方法もある。さらに、椎体・縦隔による太い帯状の“透過性低下部位”を認識し、その帯状の透過性低下部位との位置関係と、肺野の“透過性亢進部位”との位置関係から左右もしくは両肺を認識する方法を採ることも可能である。また、図２Ｈに示すように、横隔膜の下側の領域についてもこの方法を適用することができる。これにより、横隔膜の下側の部位や、心臓を認識することも可能である。
【0120】
さらに、空気は、透過性が一番高く、肺野よりも透過性が高い部位であることから、空気も考慮して計算することが望ましい。すなわち、画面上の空気の位置によって、以下のように判断することができる。
（画面の右上の空気の領域）＞（画面の左上の空気の領域）である場合は、左肺と認識する。これは、肩周りは人体外の空気の領域が撮影上広くなるからである。
（画面の左上の空気の領域）＞（画面の右上の空気の領域）である場合は、右肺と認識する。これも、上記と同様に、肩周りは人体外の空気の領域が撮影上広くなるからである。
次に、（画面の右上の空気の領域）≒（画面の左上の空気の領域）である場合は、両肺と認識する。これは、空気の領域が左右同程度であるからである。
（画面の右上の空気の領域）＞（画面の左上の空気の領域）である場合は、左肺と認識する。これは、肩周りは人体外の空気の領域が撮影上広くなるからである。
（画面の左上の空気の領域）＞（画面の右上の空気の領域）である場合は、右肺と認識する。これも、上記と同様に、肩周りは人体外の空気の領域が撮影上広くなるからである。
次に、（画面の右上の空気の領域）≒（画面の左上の空気の領域）である場合は、両肺と認識する。これは、空気の領域が左右同程度であるからである。
【0121】
なお、横隔膜下に腸管の空気が入ることがあり、その際は認識されなくなる場合がある。このため、肺野の中心部から、縦隔側、心臓側、横隔膜側など、最初に大まかな肺野およびその周りの透過性低下部位を認識し、そのラインを肺野の淵を認識することもできる。この手法は、例えば
「https://jp.mathworks.com/help/images/examples/block-processing-large-images_ja_JP.html」に開示されている技術を用いることも可能である。
「https://jp.mathworks.com/help/images/examples/block-processing-large-images_ja_JP.html」に開示されている技術を用いることも可能である。
【0122】
これにより、ある患者と他の患者との比較や数値化が可能となる。また、正常肺または正常血管と、典型的な異常肺機能または異常血流との比較や数値化が可能となる。さらに、ある患者の異なる時間における肺機能や肺血流の相対評価として、モデル化した肺および標準血流を使用することが可能となる。このようなモデル化した肺および標準血流は、種々のタイプの典型患者、健康な人の典型例を集合させ、モデル化した肺および標準血流とし、ある患者に形態的に当てはめて評価する際の指標として用いることが可能である。
【0123】
［肺野の描画］
一般的に、肺野には透過性の低い肋骨が含まれるため、「density」のみを指標として肺の輪郭を機械的に同定することは難しい。そこで、本明細書では、ベジエ曲線および直線の組み合わせを用いて肺野の輪郭を仮に描画し、合致性が高くなるように、肺輪郭を調整する手法を採用する。
一般的に、肺野には透過性の低い肋骨が含まれるため、「density」のみを指標として肺の輪郭を機械的に同定することは難しい。そこで、本明細書では、ベジエ曲線および直線の組み合わせを用いて肺野の輪郭を仮に描画し、合致性が高くなるように、肺輪郭を調整する手法を採用する。
【0124】
例えば、左肺の輪郭を４本のベジェ曲線と１本の直線で表現すると、肺輪郭上の５点と、制御点４点を求めることで、肺輪郭を描画することが可能になる。点の位置をずらして、複数の肺輪郭を描画し、“輪郭内の「density」の合計値が最大になる”、“輪郭線の内側と外側の数ピクセルの「density」合計の差分が最大になる”等の条件を用いて合致性を評価することで、肺輪郭を精度高く検出することができるようになる。実際には、比較的エッジが検出しやすい肺の上部の輪郭や、後述の方法で検出した横隔膜の位置から数点の位置を同定することも可能であり、上述のシミュレーションの試行回数を抑えることができる。古典的な二値化による輪郭抽出により、外縁に近い点を抽出し、最小二乗法等を利用して、ベジェ曲線の制御点位置を調整することも可能である。
【0125】
図３Ａおよび図３Ｂは、肺野の輪郭を、ベジエ曲線および直線の両方を用いて描画した例を示す図である。図３Ａは肺の面積が最大となる場合（極大輪郭）を示し、図３Ｂは肺の面積が最小となる場合（極小輪郭）を示す。各図において、「cp1~cp5」は制御点を示し、「p1~p5」は、ベジエ曲線上または直線上の点を示す。このように、極大輪郭と極小輪郭が把握できると、途中の輪郭を計算によって求めることが可能となる。例えば、呼気の１０％、２０％…の状態を表示することが可能となる。このように、本実施形態によれば、少なくとも一つ以上のベジエ曲線（Bezier curve）を用いて、少なくとも肺野、血管または心臓を描画することが可能となる。なお、以上の手法は、肺に限定されるわけではなく、「臓器の検出」として、他の臓器にも適用可能である。また、例えば、特定のフレームにおいて予め定められた解析範囲（腫瘍、脳の視床下部、基底核、内包の境界など）上に、少なくとも一つ以上のベジエ曲線（Bezier curve）を用いて、他のフレームにおいて解析範囲に対応する範囲を検出する処理を実行することも可能である。
【0126】
また、平面的な画像のみならず、立体的な画像（３Ｄ画像）についても適用可能である。曲面の方程式を定義し、その制御点を設定することによって、複数の曲面で囲まれた範囲を解析対象とすることが可能となる。
【0127】
［横隔膜または呼吸と連動する動的部位の動きの検出］
連続撮影された画像において、横隔膜または呼吸と連動する動的部位の動きを検出することが可能である。連続撮影された画像において、任意の間隔で画像を選択し、画像間の差分を計算すると、特にコントラストの大きい領域について、差分が大きくなる。この差分を適切に可視化することによって、動きのあった領域を検出することができる。可視化の際には、閾値によるノイズ除去や、最小二乗法などを活用したカーブフィッティング等で差分の絶対値が大きいエリアの連続性を強調することもできる。
連続撮影された画像において、横隔膜または呼吸と連動する動的部位の動きを検出することが可能である。連続撮影された画像において、任意の間隔で画像を選択し、画像間の差分を計算すると、特にコントラストの大きい領域について、差分が大きくなる。この差分を適切に可視化することによって、動きのあった領域を検出することができる。可視化の際には、閾値によるノイズ除去や、最小二乗法などを活用したカーブフィッティング等で差分の絶対値が大きいエリアの連続性を強調することもできる。
【0128】
肺野においては横隔膜や心臓が接するラインのコントラストが際立っており、図４Ａに示すように、２枚の肺画像において差分を取り、一定の閾値を設定して差分を可視化すると、図４Ｂに示すように、横隔膜や心臓が接するラインを可視化することができる。
【0129】
［横隔膜の動きの推定］
本手法では、対象画像間において、横隔膜が動いている場合は横隔膜位置を検出可能であるが、横隔膜の動きが緩やかになる箇所の検出は困難となる。すなわち、呼気吸気が切り替わるタイミングや、呼吸を止めている間、撮影の開始直後や終了直前では検出が難しい。本手法においては、任意の補完方法を用いて、横隔膜の動きを推定する。
本手法では、対象画像間において、横隔膜が動いている場合は横隔膜位置を検出可能であるが、横隔膜の動きが緩やかになる箇所の検出は困難となる。すなわち、呼気吸気が切り替わるタイミングや、呼吸を止めている間、撮影の開始直後や終了直前では検出が難しい。本手法においては、任意の補完方法を用いて、横隔膜の動きを推定する。
【0130】
前述の方法を用いて、図４Ｂに示すように横隔膜ラインを可視化した後、縦１０２４ｐｘの画像を縦８ｐｘごとに１２８個の長方形に分割し、各長方形領域に含まれる信号値を合計し、図４Ｃに示すように棒グラフ化した。複数のピークのうち、点線の矩形で示す、最も下の座標に来るピークが横隔膜の位置を示していると期待される。通常の立位ＸＰ画像において、横隔膜は曲線として表示されるが、この座標を横隔膜の位置と近似する。
【0131】
本手法で全画像について横隔膜位置を検出すると、図５に示すように「ピーク位置」が検出された。この検出した値に対して補正を行なうことで、横隔膜の動きを推定する。まず、差分が一定値よりも大きい場合は外れ値とみなして除外する（図５中の細い実線）。外れ値を除外したデータを、任意のクラスターに分割し、それぞれのクラスターに対して４次の曲線回帰を行ない、結果を繋ぎ合わせた（図５中の太い実線）。本解析では回帰分析を行なったが、本発明はこれに限定されるわけではなく、スプライン補間など任意の補完方法を用いる事が可能である。
【0132】
［動的部位検出の精緻化］
動的部位のコントラストはラインに沿って一様でない場合がある。その場合はノイズ除去に使用する閾値を変更して、複数回検出処理を行なうことによって、動的部位の形状をより正確に検出することができる。例えば、左肺において、横隔膜のラインのコントラストは、人体内部にいくに従って弱くなる傾向がある。図４Ｂにおいては、横隔膜の右半分しか検出できていない。このとき、ノイズ除去に利用した閾値の設定を変えることによって、横隔膜の左半分の残りの部分を検出することもできる。この処理を複数回繰り返すことによって、横隔膜全体の形状を検出することが可能となる。本手法によって、横隔膜の位置だけでなく、形状について線や面の変化率や変化量を数値化することも可能となり、新たな診断に役立てることができる。
動的部位のコントラストはラインに沿って一様でない場合がある。その場合はノイズ除去に使用する閾値を変更して、複数回検出処理を行なうことによって、動的部位の形状をより正確に検出することができる。例えば、左肺において、横隔膜のラインのコントラストは、人体内部にいくに従って弱くなる傾向がある。図４Ｂにおいては、横隔膜の右半分しか検出できていない。このとき、ノイズ除去に利用した閾値の設定を変えることによって、横隔膜の左半分の残りの部分を検出することもできる。この処理を複数回繰り返すことによって、横隔膜全体の形状を検出することが可能となる。本手法によって、横隔膜の位置だけでなく、形状について線や面の変化率や変化量を数値化することも可能となり、新たな診断に役立てることができる。
【0133】
このように検出された横隔膜の位置または形状を診断に利用することが可能となる。すなわち、本願発明では、横隔膜の座標をグラフ化し、上述したように計算された曲線（局面）、若しくは直線を用いて、胸郭や横隔膜の座標の計算をし、また、心拍や血管拍、肺野の「density」などを、周期に対応した位置、座標としてグラフ化したりすることが可能である。このような手法は、呼吸と連動する動的部位についても適用可能である。
【0134】
このような手法により、吸気、呼気でのHzだけでなく、横隔膜または呼吸と連動する動的部位の周波数（Hz）が変化した場合、その変化に応じた周波数帯域で計測できるようになる。そして、ＢＰＦ（band pass filter）のスペクトル抽出の際に、一定の範囲において、呼吸それぞれの状態に応じてＢＢＦを据えること、呼吸それぞれの「reconstruction phase」でＢＰＦの位置の軸が変動し、最適な状態が生じえること、それを合わせた変動性のＢＰＦを作成することが可能となる。これにより、呼吸が遅くなったり、止めたとき（Hz=0）のように、呼吸のリズムの変動があっても、それに応じた画像を提供することが可能となる。
【0135】
また、呼吸要素が呼気または吸気の全体に占める割合に基づいて、呼気または吸気の全体の周波数を計算するようにしても良い。なお、横隔膜の検出において、複数回施行し、信号または波形が安定しているものを選択するようにしても良い。以上により、検出した横隔膜の位置若しくは形状、または呼吸と連動する動的部位の位置若しくは形状、から、呼吸要素の少なくとも一つの周波数を計算することが可能となる。横隔膜または動的部位の位置または形状が把握できると、呼吸要素の周波数を把握することが可能となる。この手法によれば、波形の一部を区切ったとしても、その後の波形を追跡することができる。このため、呼吸要素の周波数が途中で変わっても、元々の呼吸要素を追従することが可能である。また、心臓の拍動などが突然変わることがあるが、心血管についても同様に適用することが可能となる。次に、本実施形態に係る各モジュールの動作について説明する。
【0136】
［呼吸機能解析］
まず、呼吸機能解析について説明する。図６Ａは、本実施形態に係る呼吸機能解析の概要を示すフローチャートである。基本モジュール１がデータベース１５からＤＩＣＯＭの画像を抽出する（ステップＳ１）。ここでは、少なくとも、一呼吸周期内に含まれる複数のフレーム画像を取得する。次に、取得した各フレーム画像において、少なくとも肺野内のある一定領域における密度（density/intensity）を用いて、呼吸要素の周期を特定する（ステップＳ２）。なお、特定した呼吸周期やこの呼吸周期から特定される波形については、以下の各ステップで用いることが可能である。
まず、呼吸機能解析について説明する。図６Ａは、本実施形態に係る呼吸機能解析の概要を示すフローチャートである。基本モジュール１がデータベース１５からＤＩＣＯＭの画像を抽出する（ステップＳ１）。ここでは、少なくとも、一呼吸周期内に含まれる複数のフレーム画像を取得する。次に、取得した各フレーム画像において、少なくとも肺野内のある一定領域における密度（density/intensity）を用いて、呼吸要素の周期を特定する（ステップＳ２）。なお、特定した呼吸周期やこの呼吸周期から特定される波形については、以下の各ステップで用いることが可能である。
【0137】
呼吸要素の周期の特定は、さらに、横隔膜の動き、胸郭の動きを用いることも可能である。また、Ｘ線の透過性が高い部位で測定される、ある一定のｖｏｌｕｍｅ、「density」／「intensity」で構成される範囲、スパイログラムなどの他の測定方法から得られるデータを用いても良い。なお、予め各臓器（ここでは肺）が有する周波数を特定しておき、その特定した周波数に対応する「density」／「intensity」を抽出しても良い。
【0138】
次に、図６Ａにおいて、肺野を自動検出する（ステップＳ３）。肺輪郭は連続的に変化するため、最大の形状と最小の形状が検出できれば、間の形状は計算により補間が可能である。ステップＳ２において特定した呼吸要素の周期に基づいて、各フレーム画像を補間することによって、各フレーム画像における肺輪郭を特定する。また、図２Ｅ～図２Ｈで示したようなパターンマッチングを行なって、肺野を検出しても良い。なお、検出した肺野について、カットオフによるノイズ除去を行なっても良い。次に、検出した肺野を複数のブロック領域に分割する（ステップＳ４）。そして、各フレーム画像における各ブロック領域の変化を計算する（ステップＳ５）。ここでは、各ブロック領域内での変化の値を平均化し、１つのデータとして表現する。
【0139】
なお、各ブロック領域内での変化の値について、カットオフによるノイズ除去を行なっても良い。次に、各ブロック領域の「density」／「intensity」の値や、また、その変化量について、上記呼吸要素の周期に基づいて、フーリエ解析または同調一致率の解析を実施する（ステップＳ６）。
【0140】
次に、フーリエ解析または同調一致率の解析により得られた結果について、ノイズ除去を行なう（ステップＳ７）。ここでは、上述したようなカットオフや、アーティファクト（artifact）の除去を行なうことができる。以上のステップＳ５からステップＳ７の動作を１回以上行ない、完了するかどうかを判断する（ステップＳ８）。ここで、ディスプレイで表示される特徴量については、合成波や他の波の混在により純度の高い要素、例えば、呼吸要素や血流要素、その他の要素の周波数同調性画像が一度のスペクトル抽出では表示できない場合もある。その際は、ディスプレイで表示される特徴量を画素値として、再度、複数回ディスプレイに至るすべてもしくは一部の解析をし直す場合がある。この作業によってさらに要素、例えば、呼吸要素や血流要素の同調性や一致性に関し純度の高い画像を取得することが可能となる。この操作については、オペレータがディスプレイの画像を視認しながら手動で行なっても良いし、出力結果からスペクトルを抽出してその分布割合を計算し直すことを自動的に行なっても良い。さらに、計算後においても、場合に応じて、ノイズカット処理、最小二乗法による穴埋め（補間）、周囲の「density」を用いた補正を行なっても良い。
【0141】
ステップＳ８において、完了しない場合は、ステップＳ５に遷移し、完了する場合は、フーリエ解析または同調一致率の解析により得られた結果を、疑似カラー画像としてディスプレイに表示する（ステップＳ９）。なお、白黒画像を表示しても良い。このように、複数のサイクルを繰り返すことによって、データの確度を高める場合もある。これにより、所望の動画を表示することが可能となる。また、ディスプレイに表示された画像を修正することで、所望の動画を得るようにしても良い。
【0142】
本実施形態では、所望の周波数または周波数帯域を計算により算出するが、実際の画像として見ると、必ずしも良い画像が表示できるとは限らない。このため、以下の手法を採る場合もある。
（１）いくつかの周波数帯を複数提示し、人的に選択する方法
（２）いくつかの周波数帯を複数提示し、ＡＩ技術によりパターン認識でよい画像を抽出する方法
（３）ＨＩＳＴＧＲＡＭの傾向、形態から選択する。すなわち、結果の信号における「Histgram」の中心部の値が高くなる傾向があり、また、動きに応じて「histgram」の値が変動するため、ＨＩＳＴＧＲＡＭの傾向、形態から選択しても良い。
（１）いくつかの周波数帯を複数提示し、人的に選択する方法
（２）いくつかの周波数帯を複数提示し、ＡＩ技術によりパターン認識でよい画像を抽出する方法
（３）ＨＩＳＴＧＲＡＭの傾向、形態から選択する。すなわち、結果の信号における「Histgram」の中心部の値が高くなる傾向があり、また、動きに応じて「histgram」の値が変動するため、ＨＩＳＴＧＲＡＭの傾向、形態から選択しても良い。
【0143】
［肺血流解析］
次に、肺血流解析について説明する。図７は、本実施形態に係る肺血流解析の概要を示すフローチャートである。基本モジュール１がデータベース１５からＤＩＣＯＭの画像を抽出する（ステップＴ１）。ここでは、少なくとも、一心拍周期内に含まれる複数のフレーム画像を取得する。次に、取得した各フレーム画像に基づいて、血管拍周期を特定する（ステップＴ２）。なお、特定した血管拍周期やこの血管拍周期から特定される波形については、以下の各ステップで用いることが可能である。血管拍周期は、上述したように、例えば、心電図や脈拍計等の他のモダリティの計測結果、心臓・肺門・主要血管など、任意の部位の「density」／「intensity」の変化を用いて血管拍を解析する。なお、予め各臓器（ここでは肺血流）が有する周波数を特定しておき、その特定した周波数に対応する「density」／「intensity」を抽出しても良い。
次に、肺血流解析について説明する。図７は、本実施形態に係る肺血流解析の概要を示すフローチャートである。基本モジュール１がデータベース１５からＤＩＣＯＭの画像を抽出する（ステップＴ１）。ここでは、少なくとも、一心拍周期内に含まれる複数のフレーム画像を取得する。次に、取得した各フレーム画像に基づいて、血管拍周期を特定する（ステップＴ２）。なお、特定した血管拍周期やこの血管拍周期から特定される波形については、以下の各ステップで用いることが可能である。血管拍周期は、上述したように、例えば、心電図や脈拍計等の他のモダリティの計測結果、心臓・肺門・主要血管など、任意の部位の「density」／「intensity」の変化を用いて血管拍を解析する。なお、予め各臓器（ここでは肺血流）が有する周波数を特定しておき、その特定した周波数に対応する「density」／「intensity」を抽出しても良い。
【0144】
次に、図７において、上述した方法で呼吸要素の周期を特定し（ステップＴ３）、その呼吸要素の周期を用いて肺野を自動検出する（ステップＴ４）。肺の輪郭の自動検出では、フレーム画像毎に、ばらつきが生ずることもあるが、ステップＴ３において特定した呼吸要素の周期に基づいて、各フレーム画像を補間することによって、各フレーム画像における肺輪郭を特定する。また、図２Ｅ～図２Ｈで示したようなパターンマッチングを行なって、肺野を検出しても良い。なお、検出した肺野について、カットオフによるノイズ除去を行なっても良い。次に、検出した肺野を複数のブロック領域に分割する（ステップＴ５）。そして、各フレーム画像における各ブロック領域の変化を計算する（ステップＴ６）。ここでは、各ブロック領域内での変化の値を平均化し、１つのデータとして表現する。なお、各ブロック領域内での変化の値について、カットオフによるノイズ除去を行なっても良い。次に、各ブロック領域の「density」／「intensity」の値や、また、その変化量について、上記血管拍周期に基づいて、フーリエ解析または同調一致率の解析を実施する（ステップＴ７）。
【0145】
次に、フーリエ解析または同調一致率の解析により得られた結果について、ノイズ除去を行なう（ステップＴ８）。ここでは、上述したようなカットオフや、アーティファクト（artifact）の除去を行なうことができる。以上のステップＴ６からステップＴ８の動作を１回以上行ない、完了するかどうかを判断する（ステップＴ９）。ここで、ディスプレイで表示される特徴量については、合成波や他の波の混在により純度の高い要素、例えば、呼吸要素や血流要素、その他の要素の周波数同調性画像が一度のスペクトル抽出では表示できない場合もある。その際は、ディスプレイで表示される特徴量を画素値として、再度、複数回ディスプレイに至るすべてもしくは一部の解析をし直す場合がある。この作業によってさらに要素、例えば、呼吸要素や血流要素の同調性や一致性に関し純度の高い画像を取得することが可能となる。この操作については、オペレータがディスプレイの画像を視認しながら手動で行なっても良いし、出力結果からスペクトルを抽出してその分布割合を計算し直すことを自動的に行なっても良い。さらに、計算後においても、場合に応じて、ノイズカット処理、最小二乗法による穴埋め（補間）、周囲の「density」を用いた補正を行なっても良い。
【0146】
ステップＴ９において、完了しない場合は、ステップＴ６に遷移し、完了する場合は、フーリエ解析または同調一致率の解析により得られた結果を、疑似カラー画像としてディスプレイに表示する（ステップＴ１０）。なお、白黒画像を表示しても良い。これにより、データの確度を高めることが可能となる。また、ディスプレイに表示された画像を修正することで、所望の動画を得るようにしても良い。
【0147】
本実施形態では、所望の周波数または周波数帯域を計算により算出するが、実際の画像として見ると、必ずしも良い画像が表示できるとは限らない。このため、以下の手法を採る場合もある。
（１）いくつかの周波数帯を複数提示し、人的に選択する方法
（２）いくつかの周波数帯を複数提示し、ＡＩ技術によりパターン認識でよい画像を抽出する方法
（３）ＨＩＳＴＧＲＡＭの傾向、形態から選択する。すなわち、結果の信号における「Histgram」の中心部の値が高くなる傾向があり、また、動きに応じて「histgram」の値が変動するため、ＨＩＳＴＧＲＡＭの傾向、形態から選択しても良い。
（１）いくつかの周波数帯を複数提示し、人的に選択する方法
（２）いくつかの周波数帯を複数提示し、ＡＩ技術によりパターン認識でよい画像を抽出する方法
（３）ＨＩＳＴＧＲＡＭの傾向、形態から選択する。すなわち、結果の信号における「Histgram」の中心部の値が高くなる傾向があり、また、動きに応じて「histgram」の値が変動するため、ＨＩＳＴＧＲＡＭの傾向、形態から選択しても良い。
【0148】
［その他の血流解析］
次に、その他の血流解析について説明する。本発明の一態様は、図１５に示すように、心臓、大動脈、肺血管、上腕動脈、頸部血管などの血流解析についても適用可能である。さらに、図示しない腹部血管や、末梢の血管などについても、同様に血流解析が可能である。図８は、本実施形態に係るその他の血流解析の概要を示すフローチャートである。基本モジュール１がデータベース１５からＤＩＣＯＭの画像を抽出する（ステップＲ１）。ここでは、少なくとも、一心拍周期内に含まれる複数のフレーム画像を取得する。次に、取得した各フレーム画像に基づいて、血管拍周期を特定する（ステップＲ２）。なお、特定した血管拍周期やこの血管拍周期から特定される波形については、以下の各ステップで用いることが可能である。血管拍周期は、上述したように、例えば、心電図や脈拍計等の他のモダリティの計測結果、心臓・肺門・主要血管など、任意の部位の「density」／「intensity」の変化を用いて血管拍を解析する。なお、予め各臓器（例えば、主要血管）が有する周波数を特定しておき、その特定した周波数に対応する「density」／「intensity」を抽出しても良い。
次に、その他の血流解析について説明する。本発明の一態様は、図１５に示すように、心臓、大動脈、肺血管、上腕動脈、頸部血管などの血流解析についても適用可能である。さらに、図示しない腹部血管や、末梢の血管などについても、同様に血流解析が可能である。図８は、本実施形態に係るその他の血流解析の概要を示すフローチャートである。基本モジュール１がデータベース１５からＤＩＣＯＭの画像を抽出する（ステップＲ１）。ここでは、少なくとも、一心拍周期内に含まれる複数のフレーム画像を取得する。次に、取得した各フレーム画像に基づいて、血管拍周期を特定する（ステップＲ２）。なお、特定した血管拍周期やこの血管拍周期から特定される波形については、以下の各ステップで用いることが可能である。血管拍周期は、上述したように、例えば、心電図や脈拍計等の他のモダリティの計測結果、心臓・肺門・主要血管など、任意の部位の「density」／「intensity」の変化を用いて血管拍を解析する。なお、予め各臓器（例えば、主要血管）が有する周波数を特定しておき、その特定した周波数に対応する「density」／「intensity」を抽出しても良い。
【0149】
次に、解析範囲を設定し（ステップＲ３）、設定した解析範囲を複数のブロック領域に分割する（ステップＲ４）。そして、各ブロック領域内での変化の値を平均化し、１つのデータとして表現する。なお、各ブロック領域内での変化の値について、カットオフによるノイズ除去を行なっても良い。次に、各ブロック領域の「density」／「intensity」の値や、また、その変化量について、上記血管拍周期に基づいて、フーリエ解析または同調一致率の解析を実施する（ステップＲ５）。
【0150】
次に、フーリエ解析または同調一致率の解析により得られた結果について、ノイズ除去を行なう（ステップＲ６）。ここでは、上述したようなカットオフや、アーティファクト（artifact）の除去を行なうことができる。以上のステップＲ５からステップＲ６の動作を１回以上行ない、完了するかどうかを判断する（ステップＲ７）。ここで、ディスプレイで表示される特徴量については、合成波や他の波の混在により純度の高い要素、例えば、呼吸要素や血流要素、その他の要素の周波数同調性画像が一度のスペクトル抽出では表示できない場合もある。その際は、ディスプレイで表示される特徴量を画素値として、再度、複数回ディスプレイに至るすべてもしくは一部の解析をし直す場合がある。この作業によってさらに要素、例えば、呼吸要素や血流要素の同調性や一致性に関し純度の高い画像を取得することが可能となる。この操作については、オペレータがディスプレイの画像を視認しながら手動で行なっても良いし、出力結果からスペクトルを抽出してその分布割合を計算し直すことを自動的に行なっても良い。さらに、計算後においても、場合に応じて、ノイズカット処理、最小二乗法による穴埋め（補間）、周囲の「density」を用いた補正を行なっても良い。
【0151】
ステップＲ７において、完了しない場合は、ステップＲ５に遷移し、完了する場合は、フーリエ解析または同調一致率の解析により得られた結果を、疑似カラー画像としてディスプレイに表示する（ステップＲ８）。なお、白黒画像を表示しても良い。これにより、データの確度を高めることが可能となる。また、ディスプレイに表示された画像を修正することで、所望の動画を得るようにしても良い。
【0152】
本実施形態では、所望の周波数または周波数帯域を計算により算出するが、実際の画像として見ると、必ずしも良い画像が表示できるとは限らない。このため、以下の手法を採る場合もある。
（１）いくつかの周波数帯を複数提示し、人的に選択する方法
（２）いくつかの周波数帯を複数提示し、ＡＩ技術によりパターン認識でよい画像を抽出する方法
（３）ＨＩＳＴＧＲＡＭの傾向、形態から選択する。すなわち、結果の信号における「Histgram」の中心部の値が高くなる傾向があり、また、動きに応じて「histgram」の値が変動するため、ＨＩＳＴＧＲＡＭの傾向、形態から選択しても良い。
（１）いくつかの周波数帯を複数提示し、人的に選択する方法
（２）いくつかの周波数帯を複数提示し、ＡＩ技術によりパターン認識でよい画像を抽出する方法
（３）ＨＩＳＴＧＲＡＭの傾向、形態から選択する。すなわち、結果の信号における「Histgram」の中心部の値が高くなる傾向があり、また、動きに応じて「histgram」の値が変動するため、ＨＩＳＴＧＲＡＭの傾向、形態から選択しても良い。
【0153】
なお、３Ｄで解析した場合、呼吸量、心拍出量、中枢の血流量を別の装置で測定することによって、相対的な値であるフーリエ解析結果から各ブロック領域の呼吸量、心拍出量、中枢の血流量を算出することが可能となる。すなわち、呼吸機能解析の場合は、呼吸量から肺換気量の推定が可能となり、肺血流解析の場合は、心（肺血管）拍出量から肺血流量の推定が可能となり、その他の血流量解析の場合は、中枢側の血流量（割合）から描出される分岐血管における推定血流量（割合）の推定が可能となる。
【0154】
また、上記のように、取得したdatabaseは、すべてを計算できればより精度の高い判断が可能であるが、コンピュータ解析を実行する上では時間を要する場合がある。そのため、任意の枚数（例えば、特定のphase）だけを抜き出し、計算を行なうようにすることもできる。これにより、解析時間を短縮することができ、さらに、呼吸の初めに観られるようなイレギュラーな個所をカットすることが可能となる。また、解析結果を表示する際に、任意の範囲を表示することができる。例えば、「呼気／吸気」の変わり目から「吸気／呼気」の変わり目の範囲を表示することによって、繰り返し再生をする際に、いわゆる「エンドレス再生」が可能となり、医師による診断をし易くすることが可能となる。
【0155】
以上説明したように、本実施形態によれば、Ｘ線動画装置で人体の画像を評価することが可能となる。デジタルデータを取得できれば、現存施設装置で概ね良好に計算可能であり、導入費用が安価で済む。例えば、Ｆｌａｔ ｐａｎｅｌ ｄｅｔｅｃｔｏｒを用いたＸ線動画装置において、被写体の検査を簡便に済ませることが可能となる。また、肺血流については、肺血栓塞栓症のスクリーニングが可能となる。例えば、Ｆｌａｔ ｐａｎｅｌ ｄｅｔｅｃｔｏｒを用いたＸ線動画装置においては、ＣＴを行なう前に本実施形態に係る診断支援プログラムを実行することによって、無駄な検査を除外することができる。また、検査が簡便であるため、緊急性の高い疾患を早期に発見し、優先的に対応することが可能となる。なお、現時点における撮影方法では、ＣＴ、ＭＲＩなどの他のmodalityでは、いくつかの課題があるが、これを解決できれば、各領域の細かい診断が可能となる。
【0156】
また、各種の血管、例えば、頸部血流狭小化のスクリーニングにも適用可能であり、また、大血管の血流評価やスクリーニングにも適用可能である。また、肺呼吸データについては、肺の部分機能検査として有効となり、肺機能検査として使用することが可能となる。また、ＣＯＰＤ、肺気腫などの疾患の同定も可能となる。さらに、術前、術後の性状の把握にも適用可能である。さらに、呼吸要素の周期および血流周期をフーリエ解析し、腹部のＸ線画像において、呼吸の波形および血流の波形を除去することで、残りの生体運動の異変、例えば、腸管イレウスなどが観察可能となる。
【0157】
なお、最初に取得した画像が、ある程度高精細である場合は、画素数が多いため、計算時間に時間がかかることもある。その場合、一定のピクセル数に画像を減らしてから計算しても良い。例えば、「4096×4096」のピクセルを実際には「1024×1024」にしてから計算することで計算時間を抑えることが可能である。
【0158】
［その他］
なお、Ｘ線画像を撮影する際に、例えば、ＡＲ法（Autoregressive Moving average model）などの予測アルゴリズムを用いることができる。呼吸要素の少なくとも一つの周波数が特定できると、この周波数に応じて、Ｘ線の照射間隔を調整するよう、Ｘ線撮影装置を制御することも可能である。例えば、呼吸要素の周波数が小さい場合（周期が長い場合）、Ｘ線撮影の回数を減らすことができる。これにより、人体の被ばく量を減らすことが可能となる。なお、頻呼吸や頻脈など、呼吸要素や心血管拍要素の周波数が大きい場合（周期が短い場合）は、照射頻度を高めて最適な画像作成を行なっても良い。
なお、Ｘ線画像を撮影する際に、例えば、ＡＲ法（Autoregressive Moving average model）などの予測アルゴリズムを用いることができる。呼吸要素の少なくとも一つの周波数が特定できると、この周波数に応じて、Ｘ線の照射間隔を調整するよう、Ｘ線撮影装置を制御することも可能である。例えば、呼吸要素の周波数が小さい場合（周期が長い場合）、Ｘ線撮影の回数を減らすことができる。これにより、人体の被ばく量を減らすことが可能となる。なお、頻呼吸や頻脈など、呼吸要素や心血管拍要素の周波数が大きい場合（周期が短い場合）は、照射頻度を高めて最適な画像作成を行なっても良い。
【0159】
また、ＤＩＣＯＭデータの保存形式であるが、圧縮をすると画像の質が低下する場合があることから、非圧縮にて保存することが望ましい。また、データの圧縮形式に応じて、計算方法を変えるようにしても良い。
【符号の説明】
【0160】
１ 基本モジュール
３ 呼吸機能解析部
５ 肺血流解析部
７ その他の血流解析部
９ フーリエ解析部
１０ 波形解析部
１１ 視覚化・数値化部
１３ 入力インタフェース
１５ データベース
１７ 出力インタフェース
１９ ディスプレイ
３ 呼吸機能解析部
５ 肺血流解析部
７ その他の血流解析部
９ フーリエ解析部
１０ 波形解析部
１１ 視覚化・数値化部
１３ 入力インタフェース
１５ データベース
１７ 出力インタフェース
１９ ディスプレイ
【図1A】
【図1B】
【図1C】
【図1D】
【図2A】
【図2B】
【図2C】
【図2D】
【図2E】
【図2F】
【図2G】
【図2H】
【図3A】
【図3B】
【図4A】
【図4B】
【図4C】
【図5】
【図6A】
【図6B】
【図6C】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】