JP5844187B2 - 画像解析装置および方法並びにプログラム - Google Patents

Description

本発明は、肺の呼吸時の動きを表す３次元画像を用いて、肺の３次元的な換気量の分布を解析する画像解析装置および方法、並びに画像解析方法をコンピュータに実行させるためのプログラムに関するものである。

近年、医療機器（例えば多検出器型ＣＴ、以下マルチスライスＣＴとする等）の進歩により質の高い３次元画像が画像診断に用いられるようになってきている。マルチスライスＣＴにおいては、Ｘ線を検出するための検出器が複数列用意されており、人体の周囲を１回転することにより、同時に複数枚の断層画像が取得可能となっている。ここで、検出器の数は技術の進歩により指数的に増大しており、現在では最大３２０列のマルチスライスＣＴも利用されてきている。また、検出器が人体の周囲を１回転するのに要する時間は、０．３５秒程度であり、さらに人体の体軸方向における撮影範囲も１６ｃｍに亘るものとなっている。このため、頭部や心臓等の部位では、複数の３次元画像を短い時間間隔により時系列で取得することが可能となっている。このように時系列で取得した３次元画像に含まれる注目する臓器を時系列順に表示することにより、心臓等の臓器が動く様子を動画像を見るように観察することが可能となる。また、頭部や心臓において、造影剤等が体内を巡る様子を、３次元的な画像において観察することも可能となる。

ところで、撮影範囲や撮影回数が多くなると被検体の被曝量も増えることになるが、近年、撮影装置等のモダリティを製造するメーカーにおいては、より低線量で高画質な画像が取得できるような画像再構成アルゴリズムの開発も行われている。このため、従来は被曝量が多くなっていた広範囲および長時間に亘る撮影も行われるようになってきている。

このような、モダリティの進歩により、肺に関しても、息を大きく吸ってから吐くという一連の呼吸の動作中において、肺領域全体（およそ３０ｃｍ程度）の３次元画像も短い時間間隔により時系列順に取得可能になりつつあり、これにより、呼吸によって肺が動く様子を動画像を見るように観察することが可能となる。

一方、呼吸器の疾患としては、気管支ぜんそくおよび肺気腫等が上げられる。従来、これらの疾患に対しては、肺活量計（スパイロメトリ）、あるいはＳＰＥＣＴの肺シンチグラム等の核医学検査を行うことにより、肺の換気量を調べることが行われている。ここで、スパイロメトリは、被験者が吐き出す息の量と吐き出す時間を計測するものであり、ＣＯＰＤ（慢性閉塞成敗疾患）等の肺疾患を診断するための検査として推奨されている。しかしながら、スパイロメトリは、肺全体の換気量や換気曲線（肺活量曲線）を求めることはできるものの、肺のどの部分の換気量が低下しているかについては確認することができない。また、ＳＰＥＣＴ検査では、肺の部位毎の換気の様子が可視化されるが、画質がそれほどよくないため、詳細な情報を得ることができない。また、放射線物質を体内に注入する必要があるため、患者の被曝量が多くなるという問題もある。

このため、呼気および吸気の２相のＣＴ画像から肺領域を抽出し、非剛体レジストレーション手法を用いて、抽出した肺領域の位置合わせを行って肺領域における変位ベクトル場を求め、その発散量を計算することにより、肺の局所的な換気量を計算する手法が提案されている（特許文献１参照）。

また、肺が動く様子を表す３次元的な時系列画像において、各３次元画像に気道を配置し、各３次元画像における気道間の３次元グリッドポイントを抽出し、抽出したポイントの体積を測定することにより、気道の軟化症の広がりを計測する手法が提案されている（特許文献２参照）。

特開２００５−２８１２１号公報 米国特許公開第２００８／０００９６８３号

しかしながら、特許文献１に記載された手法では、呼気および吸気の２相のＣＴ画像から肺の局所的な換気量を算出しているため、肺のいずれの部分から空気が取り込まれ始め、いずれの部分から空気が出ていくのかという、換気量の時間的な変化を確認することができない。一方、特許文献２に記載された手法は、呼吸時における肺の３次元的な動きおよび気道における軟化症の広がりを確認できるものの、肺の換気量を確認することはできない。

本発明は上記事情に鑑みなされたものであり、呼吸時における肺の３次元的な動きとともに、肺の局所的な異常を確認できるようにすることを目的とする。

本発明による画像解析装置は、肺についての時相が異なる一連の３次元画像に基づいて、肺の３次元的な換気量の分布を解析する画像解析装置であって、
３次元画像のそれぞれから肺領域を抽出する肺領域抽出手段と、
一連の３次元画像間において肺領域の位置合わせを行って、肺領域における変位ベクトル場を算出する位置合わせ手段と、
３次元画像のそれぞれにおいて、変位ベクトル場に基づいて、変位ベクトル場の各点における、換気量の時間的な変化を表す局所換気量関数を算出する関数算出手段と、
局所換気量関数と基準となる基準換気量関数との相違を定量化して、この相違を表す定量値を算出する定量化手段とを備えたことを特徴とするものである。

「肺についての時相が異なる一連の３次元画像」としては、同一被写体の肺を含む胸部を短い時間間隔で連続して撮影を行うことにより取得し、時系列順に表示することにより、呼吸による肺の動きを再生することが可能であれば、任意の画像を用いることができる。具体的には、ＣＴ画像あるいはＭＲ画像から生成した３次元画像等を用いることができる。

なお、本発明による画像解析装置においては、定量化手段を、局所換気量関数と基準換気量関数との差分を表す差分関数を、定量値として算出する手段としてもよい。

また、本発明による画像解析装置においては、定量化手段を、局所換気量関数と基準換気量関数との差分を表す差分関数を算出し、差分関数における換気量の絶対値の積分値または絶対値の最大値を定量値として算出する手段としてもよい。

また、本発明による画像解析装置においては、定量化手段を、局所換気量関数と基準換気量関数とにおける、所定の換気量に達するまでに時間の相違を、定量値として算出する手段としてもよい。

また、本発明による画像解析装置においては、表示される３次元画像または３次元画像から取得される２次元画像上において、定量値を可視化する表示制御手段をさらに備えるものとしてもよい。

また、本発明による画像解析装置においては、表示制御手段を、３次元画像または２次元画像のそれぞれにおける変位ベクトル場の各点に、定量値をマッピングすることにより、定量値を可視化する手段としてもよい。

また、本発明による画像解析装置においては、位置合わせ手段を、非剛体レジストレーション手法により位置合わせを行う手段としてもよい。

また、本発明による画像解析装置においては、位置合わせ手段を、３次元画像間の相関に基づいて位置合わせを行う手段としてもよい。

また、本発明による画像解析装置においては、３次元画像を、ＣＴ画像またはＭＲ画像から生成されてなるものとしてもよい。

また、本発明による画像解析装置においては、定量化手段を、局所換気量関数および基準換気量関数を、時間軸方向に正規化した後に、相違を定量化する手段としてもよい。

また、本発明による画像解析装置においては、定量化手段を、局所換気量関数および基準換気量関数を、局所換気量関数または基準換気量関数の最大値により正規化した後に、相違を定量化する手段としてもよい。

また、本発明による画像解析装置においては、基準換気量関数を、３次元画像を取得した同一被検体についての肺全体の体積変化を表す関数、同一被検体についての肺活量関数、同一被検体の肺の部分領域についての体積変化を表す関数、多数の健常者に基づいて算出した換気量関数、または数学的に算出された換気量関数としてもよい。

本発明による画像解析方法は、肺についての時相が異なる一連の３次元画像に基づいて、肺の３次元的な換気量の分布を解析する画像解析方法であって、
３次元画像のそれぞれから肺領域を抽出し、
一連の３次元画像間において肺領域の位置合わせを行って、肺領域における変位ベクトル場を算出し、
３次元画像のそれぞれにおいて、変位ベクトル場に基づいて、変位ベクトル場の各点における、時間と換気量との関係を表す局所換気量関数を算出し、
局所換気量関数と基準となる基準換気量関数との相違を定量化して、この相違を表す定量値を算出することを特徴とするものである。

なお、本発明による画像解析方法をコンピュータに実行させるためのプログラムとして提供してもよい。

本発明によれば、肺についての時相が異なる一連の３次元画像のそれぞれから肺領域が抽出され、一連の３次元画像間において肺領域の位置合わせが行われて、肺領域における変位ベクトル場が算出される。そして、３次元画像のそれぞれにおいて、変位ベクトル場に基づいて、変位ベクトル場の各点における局所換気量関数が算出され、局所換気量関数と基準となる基準換気量関数との相違が定量化されて、相違が定量値として算出される。このため、３次元画像を時系列順に表示することにより、呼吸時における肺の３次元的な動きを確認することができる。また、定量値により、肺の局所的な換気量の異常を確認することができる。

また、定量値を３次元画像上において可視化することにより、呼吸時における肺の３次元的な動きを、肺の局所的な換気量とともに確認することができる。

また、局所換気量関数および基準換気量関数を、時間軸方向に正規化した後に相違を定量化する、あるいは局所換気量関数または基準換気量関数の最大値により正規化した後に相違を定量化することにより、より正確に相違を定量化することができる。

本発明の実施形態による画像解析装置の構成を示す概略ブロック図 変位ベクトルの発散と変形による体積変化との関係を概念的に示す図 平均モデル換気量関数の算出を説明するための図 換気量関数の正規化を説明するための図 肺のある部位における正規化された局所換気量関数と正規化された基準換気量関数との比較を説明するための図 差分関数を示す図 １次元ルックアップテーブルを示す図 肺のＶＲ画像の表示を説明するための図 本実施形態において行われる処理を示すフローチャート 局所換気量関数と基準換気量関数との定量化の他の例を説明するための図 ＭＰＲ断面上の肺領域に定量値をオーバーレイ表示した状態を示す図

以下、図面を参照して本発明の実施形態について説明する。図１は本発明の実施形態による画像解析装置の構成を示す概略ブロック図である。なお、図１に示す画像解析装置１の構成は、補助記憶装置に読み込まれた画像解析プログラムをコンピュータ上で実行することにより実現される。この画像解析プログラムは、ＣＤ−ＲＯＭ等の記録媒体に記録され、もしくはインターネット等のネットワークを介して配布され、コンピュータにインストールされる。

本実施形態による画像解析装置１は、画像取得部１０、記憶部１２、肺領域抽出部１４、位置合わせ部１６、関数算出部１８、定量化部２０、表示制御部２２および入力部２４を備える。

画像取得部１０は、マルチスライスＣＴ装置またはＭＲＩ装置等のモダリティ２において、被写体の胸部を所定の時間間隔Δｔで撮影して得られた複数の３次元画像Ｖｉ（ｉ＝１〜ｎ）からなる３次元画像群Ｖｇを取得する、通信インターフェースの機能を有する。なお、３次元画像群ＶｇはＬＡＮ経由でモダリティ２から送信される。

ここで、３次元画像Ｖｉは、診断対象となる胸部を断層面に垂直な方向に沿って順に得られる２次元の断層画像を積層することによって取得されるものであり、本実施形態においては、モダリティ２において撮影された複数の断層画像を重ね合わせることにより生成される。なお、ＣＴ装置を用いて取得した３次元画像は、３次元空間上での格子点を構成するボクセル（すなわち画素位置）毎にＸ線の吸収量を蓄えたデータとなり、各画素位置に対して１つの信号値（ＣＴ装置で撮影した場合は、Ｘ線の吸収量を示す値）が与えられたデータとなる。

３次元画像群Ｖｇは、例えば、被写体を一定の時間間隔Δｔからなる異なる時相Ｔ１、Ｔ２、・・・、Ｔｎで撮影を行うことにより取得した、一連の３次元画像Ｖｉから構成される。なお、被写体の撮影は、呼気−吸気−呼気の順、または吸気−呼気−吸気の順にて行われることが好ましいが、呼気−吸気の順または吸気−呼気の順で行うようにしてもよい。

なお、３次元画像Ｖｉには、ＤＩＣＯＭ（Digital Imaging and Communications in Medicine）規格で規定された付帯情報が付加される。付帯情報は、例えば、個々の３次元画像Ｖｉにより表される３次元画像を識別するための画像ＩＤ、被写体を識別するための患者ＩＤ、検査を識別するための検査ＩＤ、画像情報毎に割り振られるユニークなＩＤ（ＵＩＤ）、その画像情報が生成された検査日、検査時刻、その画像情報を取得するための検査で使用されたモダリティの種類、患者氏名、年齢、性別等の患者情報、検査部位（撮影部位、本実施形態においては胸部）、撮影条件（造影剤の使用有無や、放射線量等）、１回の検査で複数の画像を取得したときのシリーズ番号あるいは採取番号等の情報が含まれうる。

記憶部１２は、ハードディスク等の大容量の記憶装置であり、３次元画像群Ｖｇが記憶される。なお、記憶部１２には、被写体が異なる（すなわち患者が異なる）、あるいは同一被写体で撮影時期が異なる複数の３次元画像群Ｖｇが記憶される。

肺領域抽出部１４は、３次元画像Ｖｉのそれぞれから肺領域を抽出する。肺領域を抽出する手法としては、３次元画像Ｖｉから画素毎の信号値をヒストグラム化し、肺領域をしきい値処理することにより抽出する方法、肺領域を表すシード点に基づく領域拡張法(Region Growing)等、任意の手法を用いることができる。なお、肺の境界付近に病変領域が存在する場合には、例えば特開２００８−２５３２９３号公報に記載されているように、病変領域である肺領域の凹部分を、肺領域の境界の凹部分前後における曲率を用いた曲線補間によって補間することにより、新たに肺領域に含めて抽出することが好ましい。

位置合わせ部１６は、３次元画像Ｖｉのそれぞれから抽出した肺領域における対応する画素位置を、３次元画像Ｖｉ間において位置合わせして、肺領域における変位ベクトル場を算出する。具体的には、「Boldea V, Sarrut D, Clippe S, Lung deformation estimation with non-rigid registration for radiotherapy treatment, Proceedings of MICCAI 2003, Lecture Notes in Computer Science 2003;2878:770-7」（参考文献１）に記載された手法を用いて、時相が隣接する３次元画像間において対応する画素位置を対応づけることにより位置合わせを行い、時相が隣接する２つの３次元画像間における肺領域の変位ベクトル場を算出する。参考文献１に記載された手法は、非剛体レジストレーション手法を用いて、２つのＣＴ画像から抽出した肺領域同士の位置合わせをボクセルベースで行うものであり、信号値の類似度が最大で、かつ変形による弾性エネルギが最小となる変位ベクトル場を求める手法である。なお、位置合わせの手法としては、参考文献１の手法に限定されるものではなく、時相が隣接する３次元画像間の相関を算出し、相関が最小となる変位ベクトル場を求めることにより行ってもよい。

なお、以降の説明において、位置合わせにより算出された隣接する２つの時相（Ｔｉ〜Ｔｉ＋１）間の変位ベクトル場を、２つの時相の間の時間におけるＤ（Ｔｉ）（ｘ，ｙ，ｚ）とする。ここで、（ｘ，ｙ，ｚ）は変位ベクトル場の各点の３次元座標である。また、特定の時間Ｔｋまでの変位ベクトル場の累積変化量をＦ（Ｔｋ）（ｘ，ｙ，ｚ）とする。累積変化量Ｆ（Ｔｋ）（ｘ，ｙ，ｚ）は、ｔ＝Ｔ１からｔ＝Ｔｋまで各時相間において取得される変位ベクトル場Ｄ（Ｔｉ）（ｘ，ｙ，ｚ）を加算したものとなる。

また、位置合わせを行う前に、各３次元画像Ｖｉを平滑化するようにしてもよい。具体的には、所定サイズ（例えば３×３×３）の平滑化フィルタを用いて、３次元画像Ｖｉの各画素位置の信号値の平均値を算出することにより、各３次元画像Ｖｉを平滑化すればよい。これにより、位置合わせを行う場合に、３次元画像Ｖｉに含まれるノイズの影響を低減できるため、より正確に位置合わせを行うことができる。

関数算出部１８は、位置合わせ部１６が算出した変位ベクトル場Ｄ（Ｔｉ）（ｘ，ｙ，ｚ）を用いて、単位時間当たりの変位ベクトル場の各点の容積の変位量を、局所換気量Ｐ（Ｔｉ）（ｘ，ｙ，ｚ）として算出する。具体的には、上記特許文献１に記載されたように、肺領域の変位ベクトル場の各点における発散を、肺内の局所換気量として算出する。ここで、局所換気量とは、肺内の任意の一部の領域（局所）内における換気量をいう。本実施形態において局所は、１または互いに隣接する複数のボクセルで構成される。肺領域の変位ベクトル場における各点の発散と肺内の局所換気量との関係は、以下のように説明される。

すなわち、換気は肺内における空気の非圧縮性の移動であり、空気の移動により肺実質が変形する。呼吸に伴う肺実質の変形は、肺を多孔性の弾性体とみなすことにより、弾性体の理論を適用することが可能となる。弾性体における局所変形は、局所における変位ベクトル場と変形テンソル（３×３の行列形式）を用いて表すことができ、とくに、座標軸の方向を変形テンソル場の主方向にとれば、変形テンソルは対角化され、そのトレース（対角成分の和）は、変位ベクトルの発散と一致することが知られている。また、このとき局所における体積変化率は、対角化された変形テンソルのトレースと近似することができる、すなわち、変形による体積の変化率は変位ベクトルの発散に一致するとみなすことが可能である。なお、肺血管等、空気の移動のない組織では、変形はしても体積変化は無視し得るため発散は０となる。

図２は変位ベクトル場の発散と変形による体積変化との関係を概念的に示す図である。なお、図２において、Ｄは点Ｏ１における変位ベクトル場を表し、ｄｘ，ｄｙ，ｄｚは単位立方体のｘ，ｙ，ｚ各座標軸方向の微小変位を示している。このとき、単位立方体の体積変化は、（１＋ｄｘ）・（１＋ｄｙ）・（１＋ｄｚ）−１で求まるが、これは、ｄｘ＋ｄｙ＋ｄｚで近似することができる。したがって、単位立方体の体積変化率は、変位ベクトル場の発散に一致するとみなせる。発散は任意のベクトル関数をＡ（ｘ，ｙ，ｚ）とするときに、一般にdivＡまたはハミルトンの演算子を用いて∇・Ａと表記される。ここで、本実施形態においては、ベクトル関数は変位ベクトル場Ｄ（Ｔｉ）（ｘ，ｙ，ｚ）であることから、発散divＤは、下記の式（１）により算出でき、算出した発散divＤが、変位ベクトル場の各点における、ある時相間における局所換気量Ｐ（Ｔｉ）（ｘ，ｙ，ｚ）となる。

さらに関数算出部１８は、各時相間の局所換気量Ｐ（Ｔｉ）（ｘ，ｙ，ｚ）を全ての時相について積分することにより、換気量の時間的な変化を表す局所換気量関数Ｑ（ｔ）（ｘ，ｙ，ｚ）を算出する。なお、ｔ＝Ｔ１〜Ｔｋまでの局所換気量関数Ｑ（ｔ）（ｘ，ｙ，ｚ）は、ｔ＝Ｔ１〜Ｔｋまで局所換気量Ｐ（ｔ）（ｘ，ｙ，ｚ）を加算したものとなる。

定量化部２０は、局所換気量関数Ｑ（ｔ）（ｘ，ｙ，ｚ）と基準となる基準換気量関数Ｖ（ｔ）（ｘ，ｙ，ｚ）との相違を定量化する。まず、基準換気量関数について説明する。本実施形態においては、基準換気量関数として、大域換気量関数、肺活量換気関数、部分換気量関数、平均モデル換気量関数または数学モデル換気量関数を用いる。なお、これらの関数は、あらかじめ算出されて記憶部１２に記憶されている。

大域換気量関数は、３次元画像Ｖｉを取得した同一被検体についての肺全体の体積変化を表す関数であり、変位ベクトル場における各点の局所換気量関数Ｑ（ｔ）（ｘ，ｙ，ｚ）を、肺領域内において体積積分（容積積分）することにより算出することができる。大域換気量関数は、３次元画像Ｖｉを取得した同一被検体についての肺領域全体の換気量の時間による変化を表すものとなる。

肺活量換気関数は、３次元画像Ｖｉを取得した同一被検体についてスパイロメトリを行うことにより取得することができる。肺活量換気関数は、３次元画像Ｖｉを取得した同一被検体について、実際の測定により得られた換気量の時間による変化を表すものとなる。

部分換気量関数は、３次元画像Ｖｉを取得した同一被検体についての肺野を複数の領域に分割した場合における、部分領域の体積変化を表す関数であり、変位ベクトル場における各点の局所換気量関数Ｑ（ｔ）（ｘ，ｙ，ｚ）を、部分領域毎に体積積分することにより算出することができる。ここで、ＣＯＰＤガイドラインにおいて用いられるGoddard分類においては、左右の肺を肺全体の長さに対する比率で上部、中部および下部の６つの部位に分割して診断が行われることから、Goddard分類にしたがって分割した６つ各部位のそれぞれを部分領域とすることができる。また、肺領域から葉間膜を抽出し、これにより肺を５葉に分割することにより得られる５つの部位を部分領域とすることもできる。また、肺領域から気管支、肺動脈および肺静脈を抽出し、気管支、肺動脈および肺静脈の特定の分岐毎の支配領域を算出し、その支配領域のそれぞれを部分領域とすることもできる。なお、部分領域の境界において、換気量関数が不連続であることが好ましくない場合は、部分領域毎の部分換気量関数を境界付近で補間することにより、部分領域の境界において、換気量関数を滑らかに変化させるようにしてもよい。

次に平均モデル換気量関数について説明する。図３は平均モデル換気量関数の算出を説明するための図である。まず多数の健常者について、上記と同様に３次元画像を用いて局所換気量関数を算出する。ここで、肺領域の大きさは個人により異なるため、所定の大きさの局所換気量関数空間Ａ１を設定し、各個人の肺領域を囲む領域Ｌｊ（ｊ＝１〜ｍ）を局所換気量関数空間Ａ１のサイズと一致するように正規化する。なお、局所換気量関数空間Ａ１および肺領域を囲む領域は直方体であるが、説明のために図３においては矩形にて示している。そして、多数の健常者についての変位ベクトル場における各点の局所換気量関数Ｑ（ｔ）（ｘ，ｙ，ｚ）を局所換気量関数空間Ａ１にマッピングし、各点での局所換気量関数の平均を算出する。この際、局所換気量関数空間Ａ１において、変位ベクトル場の各点が存在しない領域，すなわち、局所換気量関数空間Ａ１において肺領域が全く存在しない領域については、その領域に最も近い点における変位ベクトル場の局所換気量関数を割り当てる。そして各点での局所換気量関数の加算平均を算出し、これを平均モデル換気量関数とする。なお、各点での局所換気量関数の加算平均を，所定サイズ（例えば３×３×３）の平滑化フィルタを用いて平滑化してもよい。これにより、ノイズの影響を低減して、より正確に平均モデル換気量関数を算出できる。

また、局所換気量関数空間Ａ１への正規化に代えて、多数の健常者のうちの一人の肺領域を基準肺領域に設定し、上述した非剛体レジストレーション手法を用いて、基準肺領域に対して他の健常者の肺領域の位置合わせを行うようにしてもよい。この場合、位置合わせ後の基準肺領域に、多数の健常者についての変位ベクトル場における各点の局所換気量関数Ｑ（ｔ）（ｘ，ｙ，ｚ）をマッピングし、各点での局所換気量関数の平均を算出することにより、平均モデル換気量関数を算出すればよい。

数学モデル換気量関数は、上述した平均モデル換気量関数を、数式にて近似したものを用いる。もしくは、「“Modeling Respiratory Lung Motion - a Biophysical Approach using Finite Element Methods”, Medical Imaging 2008: Physiology, Function, and Structure from Medical Images edited by Xiaoping P. Hu, Anne V. Clough, Proc. of SPIE Vol. 6916, 69160N, (2008)」（参考文献２）、あるいは「“A 4D Model Generator of the Human Lung”, Forma, 26, 19-24, 2011」（参考文献３）に記載されたように、あらかじめシミュレーションされた肺の換気モデルから算出された肺換気量関数を使用してもよい。参考文献２に記載された手法は、モデリングされた呼吸換気量関数が、複数時相の３次元画像から抽出された特定の特徴点の位置とのずれを検証する手法である。また、参考文献３に記載された手法は、肺の人工モデルを作成し、その人工モデルから平均的な換気量関数を算出する手法である。

定量化部２０は、局所換気量関数Ｑ（ｔ）（ｘ，ｙ，ｚ）と基準換気量関数Ｖ（ｔ）（ｘ，ｙ，ｚ）との相違を定量化し、この相違を定量値として算出する。定量化を行うに際し、定量化部２０は、局所換気量関数Ｑ（ｔ）（ｘ，ｙ，ｚ）と基準換気量関数Ｖ（ｔ）（ｘ，ｙ，ｚ）とを正規化する。

ここで、基準換気量関数が、大域換気量関数および部分換気量関数のように，同一の被検体から局所換気量関数と同時に取得したものである場合、基準換気量関数と局所換気量関数との時間軸方向のずれはないため、両関数の時間軸方向の正規化は不要である。一方、定量化に際しては、両関数の形状を比較することが重要であるため、局所換気量関数Ｑ（ｔ）（ｘ，ｙ，ｚ）および基準換気量関数Ｖ（ｔ）（ｘ，ｙ，ｚ）をその全体の変量により正規化する。具体的には、局所換気量関数Ｑ（ｔ）（ｘ，ｙ，ｚ）および基準換気量関数Ｖ（ｔ）（ｘ，ｙ，ｚ）を、それぞれの最大値Ｑ（ｔ）ｍａｘおよびＶ（ｔ）ｍａｘにより正規化する。正規化した局所換気量関数Ｑ（ｔ）（ｘ，ｙ，ｚ）および基準換気量関数Ｖ（ｔ）（ｘ，ｙ，ｚ）は、それぞれＱ（ｔ）（ｘ，ｙ，ｚ）／Ｑｔ（ｍａｘ）およびＶ（ｔ）（ｘ，ｙ，ｚ）／Ｖ（ｔ）ｍａｘとなる。

一方、基準換気量関数が肺活量換気関数、平均モデル換気量関数および数学モデル換気量関数である場合、両関数の時間軸方向においても正規化が必要である。このため、量関数の呼気のタイミングおよび吸気のタイミングにより正規化を行う。図４は換気量関数の正規化を説明するための図である。なお、図４においては、基準換気量関数Ｖ（ｔ）（ｘ，ｙ，ｚ）を正規化する場合について説明する。図４に示すように、局所換気量関数Ｑ（ｔ）（ｘ，ｙ，ｚ）の吸気のタイミングＴＱ１および呼気のタイミングＴＱ２と、基準換気量関数Ｖ（ｔ）（ｘ，ｙ，ｚ）の吸気のタイミングＴＶ１および呼気のタイミングＴＱ２とはずれている。このため、基準換気量関数Ｖ（ｔ）（ｘ，ｙ，ｚ）における吸気のタイミングＴＶ１および呼気のタイミングＴＶ２が、局所換気量関数Ｑ（ｔ）（ｘ，ｙ，ｚ）の吸気のタイミングＴＱ１および呼気のタイミングＴＱ２とそれぞれ一致するように、基準換気量関数Ｖ（ｔ）（ｘ，ｙ，ｚ）を正規化する。さらに、時間軸方向に正規化した基準換気量関数を全体の変量により正規化するとともに、局所換気量関数Ｑを全体の変量により正規化する。以下、正規化した局所換気量関数Ｑ（ｔ）（ｘ，ｙ，ｚ）および基準換気量関数Ｖ（ｔ）（ｘ，ｙ，ｚ）をＱ′（ｔ）（ｘ，ｙ，ｚ），Ｖ′（ｔ）（ｘ，ｙ，ｚ）とする。

ここで、肺の各部位に正常に空気が出入りしている場合は、正規化した局所換気量関数および基準換気量関数は一致するはずである。しかしながら、実際には、特定疾患の影響、重力の影響、あるいは生理学的に正常な呼吸運動により、これらが一致しないことがある。図５は肺のある部位における正規化された局所換気量関数と正規化された基準換気量関数との比較を説明するための図である。図５において横軸は時間を、縦軸は換気量を示している。また、呼吸の開始の時間を基準（すなわち時間＝０）とし、呼気−吸気−呼気を行った状態での換気量を示している。肺の疾患の診断を行う際には、換気量が一致しない状況を確認することが重要である。このため、定量化部２０は、正規化した局所換気量関数および基準換気量関数の差分値の関数である差分関数Ｂ（ｔ）（ｘ，ｙ，ｚ）（＝Ｑ′（ｔ）（ｘ，ｙ，ｚ）−Ｖ′（ｔ）（ｘ，ｙ，ｚ））を、定量値として算出する。図６は差分関数を示す図である。図６に示すように、差分関数Ｂ（ｔ）は、時間の経過による局所換気量と基準換気量との相違、より具体的には、局所位置の換気が基準となる換気よりもどの程度遅れているかを表すものとなる。

なお、基準換気量関数が部分換気量関数である場合には、その部分換気量関数を取得した部分領域についての局所換気量関数とその部分換気量関数とを用いて差分関数を算出すればよい。

表示制御部２２は、３次元画像Ｖｉをディスプレイ４に表示するに際し、差分関数Ｂ（ｔ）（ｘ，ｙ，ｚ）を３次元画像Ｖｉとともに表示する。具体的には、３次元画像Ｖｉから肺領域を抽出して、肺のボリュームレンダリング画像またはサーフェスレンダリング画像（以下、単にＶＲ画像とする）を表示する。さらに、差分関数Ｂ（ｔ）（ｘ，ｙ，ｚ）を、肺のＶＲ画像の表面にマッピング表示する。すなわち、各時相における肺の表面の画素位置における信号値を、その画素位置に対応する差分関数Ｂ（ｔ）（ｘ，ｙ，ｚ）の値に応じた色に変換して表示する。この色の変換は、例えば図７に示すように、横軸に差分関数Ｂ（ｔ）の値を、縦軸に色（Ｒ，Ｇ，Ｂ）が設定された１次元のルックアップテーブルを用いて行う。なお、図７においては１つのルックアップテーブルのみを示しているが、実際には、Ｒ，Ｇ，Ｂの各色について３つのルックアップテーブルが用意されている。

表示制御部２２は、ルックアップテーブルを参照して、肺の表面の各画素の信号値を、Ｒ，Ｇ，Ｂの各色からなる表示画素値に変換する。そして、変換された肺のＶＲ画像をディスプレイ４に時系列順に表示する。図８は肺のＶＲ画像の表示を説明するための図である。図８に示すように、ディスプレイ４には、肺における換気量が異常となる部分が、肺の呼吸に応じた変形の状態とともに表示される。なお、ディスプレイ４には肺の３次元画像が表示されるが、図８においては説明のために２次元画像にて示している。また、図８においては、吸気の状態において、局所的に換気量の異常が見られる部位を斜線にて示している。

入力部２４は、キーボードおよびマウス等の公知の入力装置からなる。

次いで、本実施形態において行われる処理について説明する。図９は本実施形態において行われる処理を示すフローチャートである。なお、複数の肺の３次元画像Ｖｉは、画像取得部１０により取得されて、記憶部１２に記憶されているものとする。操作者が入力部２４を操作することにより、表示する３次元画像の選択がなされると（ステップＳＴ１肯定）、肺領域抽出部１４が、選択された３次元画像Ｖｉを記憶部１２から読み出し、３次元画像Ｖｉのそれぞれから肺領域を抽出する（ステップＳＴ２）。そして、位置合わせ部１６が、３次元画像Ｖｉのそれぞれから抽出した肺領域に対応する画素位置を、３次元画像Ｖｉ間において位置合わせする（ステップＳＴ３）。これにより、３次元画像Ｖｉの各時相における変位ベクトル場が算出される。

そして、関数算出部１８が、変位ベクトル場の各点における換気量の時間的な変化を表す局所換気量関数Ｑ（ｔ）（ｘ，ｙ，ｚ）を算出する（ステップＳＴ４）。次いで、定量化部２０が定量化を行って、局所換気量関数および基準換気量関数の差分値の関数である差分関数Ｂ（ｔ）（ｘ，ｙ，ｚ）を、局所換気量関数および基準換気量関数の相違を表す定量値として算出する（定量化、ステップＳＴ５）。さらに、表示制御部２２が、差分関数がマッピングされた肺のＶＲ画像をディスプレイ４に時系列順に表示し（ステップＳＴ６）、処理を終了する。

ここで、癌、気管支ぜんそくおよび肺梗塞等の肺疾患では、肺の一部が硬化してしまい、肺胞の動きがよくない部分が存在する。しかしながら、単に画像を見るだけではそのような部分の診断は難しい。このため、スパイロメトリのように換気量を算出する検査方法もあるが、換気量だけでは、最終的に肺がどの程度空気を取り込んだかが分かるのみであり、肺の動きの状態は分からない。

本実施形態においては、局所換気量関数と基準換気量関数との相違を表す差分関数を３次元画像にマッピングして時系列順に表示するようにしたため、呼吸時における肺の３次元的な動きを確認することができ、かつマッピングされた差分関数により、周囲の状態と比較しての肺の局所的な異常を確認することができる。なお、肺野の一部が完全に硬化してしまって、空気の換気が行われていない場合、換気量が原理的には０になる。本実施形態においては、所定値以上の換気量を有する部位において、換気量が周囲の部位の換気量と比較して正常か異常かを確認でき、また、異常の場合は、どの程度異常かを目で見て確認することができる。したがって、本実施形態を用いて肺の状態を定期的に経過観察することにより、病状が悪化傾向にあるかどうかの診断を容易に行うことができる。

ここで、基準換気量関数を大域換気量関数とした場合、仮に肺全体に空気が均等に取り込まれて呼吸が行われたとすると、肺の各部位における局所換気量関数は、基準換気量関数と類似する関数形状となるはずである。しかしながら、肺に何らかの疾患がある場合は、疾患のある部位の局所換気量関数は、基準換気量関数と異なる関数形状となる。このため、基準換気量関数を大域換気量関数とすることにより、肺の局所的な異常を検出することができる。

一方、基準換気量関数を大域換気量関数とした場合、サンプリング点、すなわち変位ベクトル場を算出した点が少ないと、正確な換気量を推定できないおそれがある。一方、スパイロメトリにおいては、容易に連続した時間の換気量関数を取得することができるため、基準換気量関数を正確かつ容易に算出できる。

また、肺領域の全体から１つの基準換気量関数を算出し、それと局所換気量関数との相違を定量値とすると、重力の影響による部位毎の換気量関数との相違、あるいは本来生理学的に正常だと判断されるべき部位毎の換気量関数との相違を区別することができず、その結果、このような相違を有する部位も異常部位として検出されてしまうおそれがある。このため、基準換気量関数を部位換気量関数とし、部位毎に部位換気量関数と局所換気量関数との相違を定量化することにより、局所的に異常のある部位のみを検出することができる。

また、基準換気量関数を、健常者を用いて事前に計算された平均モデル換気量関数および数学モデル換気量関数とすることにより、さらに異常部位の特定が容易となるとともに、肺野の重力による換気量の偏り、あるいは生理学的に存在する換気量の偏りの影響を排除した異常部位の特定が可能となる。

なお、上記実施形態においては、肺のＶＲ画像における肺の表面に差分関数Ｂ（ｔ）をマッピングしているが、ＶＲ画像における肺を切断した断面上に差分関数Ｂ（ｔ）をマッピングしてもよい。

また、上記実施形態においては、差分関数Ｂ（ｔ）を定量値として算出しているが、図６に示す差分関数Ｂ（ｔ）の絶対値の最大値｜Ｂ（ｔｋ）｜を定量値として算出してもよい。この絶対値の最大値｜Ｂ（ｔｋ）｜は、その差分関数Ｂ（ｔ）を算出した部位における、正常な部位との最大の換気量の相違を表すものとなる。また、絶対値｜Ｂ（ｔ）｜が最大となる時間Ｔｋを定量値としてもよい。この時間Ｔｋは、その差分関数Ｂ（ｔ）を算出した部位における正常な部位との換気量の相違が最大となる時間を表すものとなる。なお、時間は例えば呼気の開始時のように所定の時点を基準とすればよい。また、差分関数Ｂ（ｔ）の積分値を定量値としてもよい。この場合、定量値は、その差分関数Ｂ（ｔ）を算出した部位における正常な部位との換気量の相違の総和を表すものとなる。

また、図１０に示すように、空気の換気が他の部位と比較して早く行われている部位については、局所換気量関数が所定の換気量Ｖｋとなるまでの時間Ｔｐ１は、基準換気量関数が所定換気量となるまでの時間Ｔｂよりも早くなる。逆に、空気の換気が他の部位と比較して遅く行われている部位については、局所換気量関数が所定の換気量となるまでの時間Ｔｐ２は、基準換気量関数が所定換気量となるまでの時間Ｔｂよりも遅くなる。このため、局所換気量関数において、ある基準となる時間（例えば図１０に示すように呼気の開始時の時間）から所定の換気量Ｖｋとなるまでの時間と、基準換気量関数における所定の換気量Ｖｋとなるまでの時間との差分値（例えばＴｐ−Ｔｂ）を定量値として算出してもよい。この場合、ある部位における定量値が負となれば、その部位は空気の換気が他の部位と比較して早く行われており、定量値が正となれば、空気の換気が他の部位と比較して遅く行われていることが分かる。

なお、このような差分関数でない定量値を算出した場合においても、肺のＶＲ画像に定量値をマッピングすればよい。この場合、時間の経過によっても肺の表面の色は変化しないこととなるが、肺の動きの状態と換気量が異常値となる部位とを確認することができる。また、この場合においても、ＶＲ画像における肺を切断した断面上に定量値をマッピングしてもよい。さらに、図１１に示すように、３次元画像ＶｉにおけるＭＰＲ断面上の肺領域に定量値をオーバーレイ表示してもよい。なお、図１１においてはオーバーレイ表示した部分を斜線にて示している。

１ 画像解析装置
２ モダリティ
４ ディスプレイ
１０ 画像取得部
１２ 記憶部
１４ 肺領域抽出部
１６ 位置合わせ部
１８ 関数算出部
２０ 定量化部
２２ 表示制御部
２４ 入力部

Claims (13)

1. 肺についての時相が異なる一連の３次元画像に基づいて、前記肺の３次元的な換気量の分布を解析する画像解析装置であって、
   前記３次元画像のそれぞれから肺領域を抽出する肺領域抽出手段と、
   前記一連の３次元画像間において前記肺領域の位置合わせを行って、前記肺領域における変位ベクトル場を算出する位置合わせ手段と、
   前記３次元画像のそれぞれにおいて、前記変位ベクトル場に基づいて、該変位ベクトル場の各点における、換気量の時間的な変化を表す局所換気量関数を算出する関数算出手段と、
   前記局所換気量関数と基準となる基準換気量関数との相違を定量化して、該相違を表す定量値を算出する定量化手段と、
   時系列順に表示される前記３次元画像または該３次元画像から取得されて時系列順に表示される２次元画像上において、前記定量値を可視化する表示制御手段とを備えたことを特徴とする画像解析装置。
2. 前記定量化手段は、前記局所換気量関数と前記基準換気量関数との差分を表す差分関数を前記定量値として算出する手段であることを特徴とする請求項１記載の画像解析装置。
3. 前記定量化手段は、前記局所換気量関数と前記基準換気量関数との差分を表す差分関数を算出し、該差分関数における換気量の絶対値の積分値または該絶対値の最大値を前記定量値として算出する手段であることを特徴とする請求項１記載の画像解析装置。
4. 前記定量化手段は、前記局所換気量関数と前記基準換気量関数とにおける、所定の換気量に達するまでの時間の相違を前記定量値として算出する手段であることを特徴とする請求項１記載の画像解析装置。
5. 前記表示制御手段は、時系列順に表示される前記３次元画像または前記２次元画像のそれぞれにおける前記変位ベクトル場の各点に、前記定量値をマッピングすることにより、前記定量値を可視化する手段であることを特徴とする請求項１から４のいずれか１項記載の画像解析装置。
6. 前記位置合わせ手段は、非剛体レジストレーション手法により前記位置合わせを行う手段であることを特徴とする請求項１から５のいずれか１項記載の画像解析装置。
7. 前記位置合わせ手段は、前記３次元画像間の相関に基づいて前記位置合わせを行う手段であることを特徴とする請求項１から５のいずれか１項記載の画像解析装置。
8. 前記３次元画像は、ＣＴ画像またはＭＲ画像から生成されてなることを特徴とする請求項１から７のいずれか１項記載の画像解析装置。
9. 前記定量化手段は、前記局所換気量関数および前記基準換気量関数を、時間軸方向に正規化した後に、前記相違を定量化する手段であることを特徴とする請求項１から８のいずれか１項記載の画像解析装置。
10. 前記定量化手段は、前記局所換気量関数および前記基準換気量関数を、前記局所換気量関数または前記基準換気量関数の最大値により正規化した後に、前記相違を定量化する手段であることを特徴とする請求項１から９のいずれか１項記載の画像解析装置。
11. 前記基準換気量関数は、前記３次元画像を取得した同一被検体についての肺全体の体積変化を表す関数、前記同一被検体についての肺活量関数、前記同一被検体の前記肺の部分領域についての体積変化を表す関数、多数の健常者に基づいて算出した換気量関数、または数学的に算出された換気量関数であることを特徴とする請求項１から１０のいずれか１項記載の画像解析装置。
12. 肺についての時相が異なる一連の３次元画像に基づいて、前記肺の３次元的な換気量の分布を解析する画像解析方法であって、
    前記３次元画像のそれぞれから肺領域を抽出し、
    前記一連の３次元画像間において前記肺領域の位置合わせを行って、前記肺領域における変位ベクトル場を算出し、
    前記３次元画像のそれぞれにおいて、前記変位ベクトル場に基づいて、該変位ベクトル場の各点における、時間と換気量との関係を表す局所換気量関数を算出し、
    前記局所換気量関数と基準となる基準換気量関数との相違を定量化して、該相違を表す定量値を算出し、
    時系列順に表示される前記３次元画像または該３次元画像から取得されて時系列順に表示される２次元画像上において、前記定量値を可視化することを特徴とする画像解析方法。
13. 肺についての時相が異なる一連の３次元画像に基づいて、前記肺の３次元的な換気量の分布を解析する画像解析方法をコンピュータに実行させるためのプログラムであって、
    前記３次元画像のそれぞれから肺領域を抽出する手順と、
    前記一連の３次元画像間において前記肺領域の位置合わせを行って、前記肺領域における変位ベクトル場を算出する手順と、
    前記３次元画像のそれぞれにおいて、前記変位ベクトル場に基づいて、該変位ベクトル場の各点における、時間と換気量との関係を表す局所換気量関数を算出する手順と、
    前記局所換気量関数と基準となる基準換気量関数との相違を定量化して、該相違を表す定量値を算出する手順と、
    時系列順に表示される前記３次元画像または該３次元画像から取得されて時系列順に表示される２次元画像上において、前記定量値を可視化する手順とをコンピュータに実行させることを特徴とするプログラム。