JP6402422B2 - 医用画像処理装置及び医用画像処理方法 - Google Patents

Description

本発明の実施形態は、医用画像処理装置及び医用画像処理方法に関する。

慢性閉塞性肺疾患（COPD：Chronic Obstructive Pulmonary Disease）の診断に関するもので、特に胸部CT画像を用いた画像診断装置に使用されるものである。COPDは慢性気管支炎と肺気腫を包含する疾患で、特に米国においては６５〜８４歳において３大疾患（がん、心疾患、脳卒中）に次ぐ第４の死因であり、２０２０年までには脳卒中を抜いて第３の死因になると推定されている。慢性気管支炎は気管や気管支が慢性的に炎症を起こし、気管支内における持続性あるいは反復性の粘液分泌の過剰状態を指し、気管支構造としては線維組織により内腔が狭められてしまう。また、肺気腫は末梢の肺胞・肺胞道が破壊されて融合し、気腔（ブラ:bulla）が生じて弾力性を失った状態となる。両方ともに、肺胞の張力で広がっていた末梢の細気管支が閉塞し、呼気能が低下する疾患である。

COPDの臨床診断は呼吸機能検査装置であるスパイロメータで測定する呼気流量を用いたスパイロメトリーによりで行われている。これは気管支拡張薬投与後において努力性肺活量（FVC）に対する１秒後の呼気量（１秒量：FEV1）の比（１秒率：ＦＥＶ１％）が７０％未満であればCOPDと診断し、更に標準的なFEV1に対するFEV1の比（％１秒率 ：％ＦＥＶ１）よりその重症度を診断している。重症度はGOLD（Global Initiative for Chronic Obstructive Lung Disease）と呼ばれる国際プロジェクトで定められたGOLD値（GOLD 0(リスク群)からGOLD４（最重症））により表される。このスパイロメトリーによるCOPD診断は肺野全体の機能を診ることに留まるため、治療戦略や治療効果判定・フォローアップのためには医用画像により画素単位毎、あるいは肺葉や肺区域単位などの肺領域毎の局所的なあるいは局在的なCOPD診断が望まれている。

そこで、吸気時の胸部ＣＴ画像を用いてCOPDによる上記の気管支構造変化やブラの存在を低濃度吸収領域（LAA：Low Attenuation Area）を提示することで局所的なCOPDの程度を表すことが既に実現されている。しかしながら、この手法ではあるＣＴ値におけるある閾値（例えば-950HU）での判定であることから、ＣＴ値そのものが常にどのＣＴ装置でも同一とは必ずしもならいないため正確な診断とは言えない。この問題に対して最近、観察する方法として胸部ＣＴ画像の呼気時画像に吸気時画像の位置合わせを行い、その位置合わせ情報から得られるボクセル毎の移動量を提示する方法、あるいは位置合わせされたボクセル毎の呼気時と吸気時のＣＴ値から判断する方法（非特許文献１）が提案されている。後者の方法の場合は呼気時のＣＴ値（ｘ軸）と吸気時のＣＴ値（y軸）を２軸とする2次元座標空間において吸気時画像に対して-950HU、呼気時画像に対して-856HUを閾値と定めて、形成される４つの象限に対して正常、COPD(実質的には肺気腫)、それらの中間的な疾患状態である機能的小気管支疾患（fSAD）を定めて画素単位で判定するものである。

しかしながら、上記の胸部ＣＴ画像の呼気時と吸気時の画像を利用したボクセルの動きやＣＴ値だけではCOPDの程度を正確に判断することは困難である。

C. J Galban,et.al.," Computed tomography-based biomarker provides unique signature for diagnosis of COPD phenotypes and disease progression", Nature Med., Vol.18, No.11, 1711-1716, 2013

目的は、ＣＴやＭＲＩによる胸部画像を用いて慢性閉塞性肺疾患診断能の向上に有効な支援情報を提供することにある。

本実施形態に係る及び医用画像処理装置は、被検体胸部の形態を表す呼吸時相の異なる複数の画像のデータを扱う。複数の画像間での部位の動き量が画素又は領域ごとに計算される。動き量と、画像の画素値から求められる特徴値と、領域の大きさの変化率との少なくとも２つの組み合わせにより慢性閉塞性肺疾患の重症度に関するレベルが画素又は領域ごとに決定される。決定されたレベルに関する情報が出力される。

図１は本実施形態に係る医用画像処理装置をＸ線コンピュータ断層撮影装置とともに示す図である。 図２は本実施形態に係る医用画像処理の手順を示す流れ図である。 図３は図２の工程Ｓ１にかかる画像例を示す図である。 図４は図２の工程Ｓ２にかかる画像例を示す図である。 図５は図２の工程Ｓ４の区域例を示す図である。 図６は図２の工程Ｖ３，Ｒ３によるレベル特定処理の変形例を示す図である。 図７は図２の工程Ｖ３，Ｒ３によるレベル特定処理の変形例を示す図である。 図８は図２の工程Ｖ３，Ｒ３によるレベル特定処理の変形例を示す図である。 図９は図２の工程Ｖ３，Ｒ３によるレベル特定処理の変形例を示す図である。 図１０は図２の工程Ｖ３，Ｒ３によるレベル特定処理の一例を示す図である。 図１１は図２の工程Ｖ３，Ｒ３によるレベル特定処理の他の例を示す図である。 図１２は図２の工程Ｖ３，Ｒ３によるレベル特定処理の他の例を示す図である。 図１３は図２の工程Ｖ３，Ｒ３によるレベル特定処理の他の例を示す図である。 図１４は図２の工程Ｖ４にかかる画像例を示す図である。 図１５は図２の工程Ｒ４にかかる画像例を示す図である。 図１６は本実施形態において重症度レベルの決定に用いるパラメータの一つである体積変化率を示す図である。 図１７は図１６の体積変化率の詳細を示す図である。 図１８は図１６の体積変化率を用いた重症度レベルの特定処理の一例を示す図である。 図１９は図１６の体積変化率を用いた重症度レベルの特定処理の他の例を示す図である。 図２０は本実施形態において重症度レベルの決定に用いるパラメータの一つである３時相の動き量を示す図である。

本実施形態に係る及び医用画像処理装置は、被検体胸部の形態を表す呼吸時相の異なる複数の画像のデータを扱う。複数の画像間での部位の動き量が画素又は領域ごとに計算される。動き量と、画像の画素値から求められる特徴値と、領域の大きさの変化率との少なくとも２つの組み合わせにより慢性閉塞性肺疾患の重症度に関するレベルが画素又は領域ごとに決定される。決定されたレベルに関する情報が出力される。

以下、図面を参照しながら本実施形態に係る医用画像処理装置を説明する。以下説明においては、被検体胸部内部の３次元形態（３次元構造）を示すボリュームデータを例示する。当該ボリュームデータは、Ｘ線コンピュータ断層撮影装置又は磁気共鳴イメージング装置により発生される。なお、本実施形態に係る医用画像処理装置は、被検体胸部内部の２次元形態（２次元構造）を示す断層像を処理対象として用いてもよい。以下単に「画像」という場合には「ボリュームデータ又は２次元画像（断層像）のデータ」を意味するものとする。また単に「画素」という場合には「ボリュームデータを構成するボクセル、又は２次元画像（断層像）を構成する画素を意味するものとする。

図１に示すように、本実施形態に係る医用画像処理装置は、それ単独で、又はＸ線コンピュータ断層撮影装置に組み込まれる。Ｘ線管１０１とＸ線検出器１０３とが一体として被検体の周囲を回転する回転／回転方式として説明するが、リング上に多数のＸ線検出器１０３が配置され、Ｘ線管１０１のみが被検体の周囲を回転する固定／回転方式と、複数のＸ線管１０１がリング上に配置され、複数のＸ線検出器１０３も同様にリング上に配置された固定／固定方式等様々な方式であってもよい。また回転／回転方式に関して、一対のＸ線管１０１とＸ線検出器１０３とが回転フレームに搭載された一管球型として説明するが、Ｘ線管１０１とＸ線検出器１０３との対が回転フレームに複数搭載されたいわゆる多管球型であってもよい。Ｘ線検出器１０３に関して、被検体を透過したＸ線をシンチレータ等の蛍光体で光に変換してからフォトダイオード等の光電変換素子で電荷に変換する間接変換形と、Ｘ線による半導体内の電子正孔対の生成及びその電極への移動すなわち光導電現象を利用した直接変換形とがあるが、いずれの方式を採用してもよい。

架台部（ガントリともいう）１００は、回転リング１０２を有する。架台駆動部１０７はホストコントローラ１１０の制御のもとで回転リング１０２を回転する。回転リング１０２にはＸ線管１０１とＸ線検出器１０３とが搭載されている。Ｘ線検出器１０３は、被検体が天板スライドにより挿入される撮影領域Ｓを挟んで、Ｘ線管１０１に対向する。高電圧発生部１０９はＸ線管１０１の両極間に印加する管電圧を発生し、陰極フィラメントに供給する管電流を発生する。高電圧発生部１０９からスリップリング１０８を介してＸ線管１０１へ管電圧が印加され且つ管電流が供給されるときＸ線管１０１からＸ線が発生される。Ｘ線管１０１で発生されたＸ線は図示しないコリメータ（Ｘ線絞り器ともいう）で例えばコーンビーム形に成形される。Ｘ線検出器１０３はコーンビーム形Ｘ線を受光するために多列型検出器（マルチスライス型検出器ともいう）が採用される。Ｘ線検出器１０３には、一般的にＤＡＳ(data acquisition system) と呼ばれているデータ収集回路１０４が接続されている。データ収集回路１０４には、Ｘ線検出器１０３の各チャンネルの電流信号を電圧に変換するＩ−Ｖ変換器と、この電圧信号をＸ線の曝射周期に同期して周期的に積分する積分器と、この積分器の出力信号を増幅するアンプと、このプリアンプの出力信号をディジタル信号に変換するアナログ・ディジタル・コンバータとが、チャンネルごとに設けられている。

データ収集回路１０４には、非接触データ伝送装置１０５を介して前処理部１０６が接続される。前処理部１０６は、データ収集回路１０４から受信したデータに対して、チャンネル間の感度不均一を補正したり、またＸ線強吸収体、主に金属部による極端な信号強度の低下又は信号脱落を補正する等の前処理を実行する。記憶装置１１２は、前処理部１０６で前処理を受けた投影データ、その他、再構成されたボリュームデータや断層像データ等を記憶する。再構成処理部１１８は、記憶された投影データに基づいてボリュームデータや断層像データを再構成する。表示部１２２は、主に画像又は後述するカラー画像を表示する。

本実施形態に係る医用画像処理装置は、呼吸時相の異なる複数の画像を処理する。典型的には、肺野が最小容積を示す呼気時相で収集された第１のボリュームデータと肺野が最大容積を示す吸気時相で収集された第２のボリュームデータとが適用される。

画像位置合わせ処理部１１９は、第１のボリュームデータと第２のボリュームデータとを解剖学上の同じ対象物どうしを同じ位置に整合させ（位置合わせ）、実際には呼気時の第１のボリュームデータと、吸気時の第２のボリュームデータとを位置合わせ（レジストレーション）する。更にこの画像位置合わせ処理部１１９は、呼気時ボリュームデータ又は吸気時ボリュームデータと標準的な呼気時又は吸気時の胸部ボリュームデータ間で解剖学上の位置を合わせる位置合わせ処理にも使用する。この位置合わせ処理は、いわゆる非剛体位置合わせ（non-rigid registration）の技術を用いてボクセルレベルの位置合わせを行う。当該非剛体位置合わせとは、被検体の呼吸による肺野の変位及び変形を考慮して解剖学的に同じ画素（ボクセル）どうしを対応付けることを言う。つまり、非剛体位置合わせ処理として、ボリュームデータ上の対象物どうしの位置合わせを行うとともに、対象物どうしのひずみを補正する処理と言える。実際的には肺領域抽出、気管支領域抽出、あるいは肺動脈の領域抽出を行い、それらの領域の呼気時、吸気時画像の差から呼気時画像空間における領域毎の変位方向と距離とを示すワープフィールド（warp field）を作成して、呼気時画像を変形する手法などを用いる。基本的には特開２０１２−１７９３５９号公報に開示された非剛体位置合わせ技術が適用される。

領域分割処理部１１７は、例えば含気肺領域に対応するＣＴ値を閾値として膨張収縮処理などにより第１、第２のボリュームデータ各々から肺野領域を抽出し、抽出した肺野領域を複数の領域に分割する。例えば肺野領域は左肺領域と右肺領域とに分割される。左右肺領域それぞれはさらに肺葉間膜に対応するＣＴ値を閾値として右肺は上葉領域、中葉領域、下葉領域の３つの肺葉領域、左肺は上葉領域、下葉領域の２つの肺葉領域に分割される。肺葉領域それぞれは気管支や気管支に伴走する肺動脈の分枝構造を利用したパターンマッチング技術等の処理により同定される葉気管支の分岐により、さらに小さな単位である区域領域に分けられる。典型的には右肺は１０個、左肺は９個の区域領域に分割される。なお、後述のベクトル処理等の対象として左右肺領域、肺葉領域、区域領域、ボクセル単位、又は近傍する所定数のボクセルからなる集合体（例えば立方体など）としてのブロック領域いずれを選択するかは入力装置１１５を介した操作者等の任意である。

ベクトル計算処理部１２１では呼吸期の異なる第１、第２のボリュームデータ間でボクセル毎に対象物の変位方向と変位距離、つまり動きベクトルを個々に計算する。処理対象が上記ブロック領域又は肺葉領域等の分割領域であるとき、ベクトル計算処理部１２１は、各領域内に存する複数のボクセルに係る複数の動きベクトルから単一の動きベクトルを発生する。例えば、各領域の複数の動きベクトルを対象とした平均ベクトルが求められ、または複数の動きベクトルの中から変位距離が中間値又は中央値を示す代表的な動きベクトルが選択される。

領域評価値計算部１２０は、上記ブロック領域毎又は肺葉領域等の分割領域毎に、各領域のＣＴ値分布の特徴を示す評価値を計算する。評価値（特徴値）として典型的には、各領域内のボクセル総数に対する空気やガスに対応する−９５０ＨＵ未満を示すボクセル数の割合（低吸収領域率（%LAV:Low Attenuation Volume）という）、各領域のＣＴ値分布（ＣＴ値ヒストグラム）から得られるＣＴ値の平均値、分散値、歪度（skewness）、尖度（kurtosis）などが上げられる。他の例として特徴値としては、画素値、画素値の差、領域内の画素数に対する閾値未満の画素値を有する画素の数の割合、領域内の画素値分布の分散、または領域内の画素値の平均値が用いられる。いずれの評価値を処理対象として選択するかは操作者等の任意である。

領域評価値計算部１２０は、さらに第１、第２のボリュームデータ間での領域の大きさの変化率、３次元では体積変化率、２次元では面積変化率を計算する。図１６に示すように、領域Ａの体積変化率ＶＲAは、例えば
|Ｖ１−Ｖ２|／Ｖ１
で与えられる。図１７に示すように、Ｖ１は、第１のボリュームデータにおける近傍する複数の画素（ベクトル）で区画される領域Ａの体積（又は面積）を示す。領域Ａは、３次元であれば、例えば近傍２７画素の範囲に設定される。なお、領域Ａは、２次元であれば、近傍９画素の範囲に設定される。この２７画素には個々に動きベクトルが求められている。第１のボリュームデータ上の２７画素は、動きベクトルを介して第２のボリュームデータ上の２７画素にそれぞれ対応する。呼気時の領域Ａは、吸気時には拡大する。第２のボリュームデータ上の２７画素は、離散している。Ｖ２は、第２のボリュームデータ上の２７画素で囲まれる領域Ａの体積（又は面積）を示す。

重症度レベル決定部１２３は、いずれかの特徴値、ボクセル（画素）又は領域ごとの動きベクトル、領域の大きさの変化率との少なくとも２つの組み合わせにより慢性閉塞性肺疾患の重症度に関するレベルを画素又は領域ごとに決定する。具体的には後述する。

カラー画像生成処理部１２４は、重症度レベル決定部１２３でボクセル毎又は４×４×４のブロック領域に又は分割領域毎に決定された重症度レベルの３次元分布を、重症度レベルに応じて予め割り当てられているカラーに係る第１又は第２のボリュームデータの解剖学的位置に一致する３次元カラー画像を生成する。さらにカラー画像生成処理部１２４は、上記３次元カラー画像から、ＭＰＲ（断面変換処理）又はボリュームレンダリング処理等により、表示部１２２で表示可能な２次元のカラー画像に変換する。

図２には、本実施形態に係る医用画像処理手順を示している。まず、図３に示すように、吸気時（吸気息止め時）に収集した投影データに基づいて３次元の胸部ＣＴ画像（第２ボリュームデータともいう）を再構成（撮像）する（Ｓ１）。肺野では成人の場合、体軸方向に４０ｃｍの範囲で撮像する必要があり、一般にはマルチスライス対応の２次元検出器１０３によりＭSＣＴスキャニングとヘリカルスキャンとを併用する。同様に、図４に示すように、呼気時（呼気息止め時）に収集した投影データに基づいて３次元の胸部ＣＴ画像（第１ボリュームデータ）を再構成する（Ｓ２）。

次に、吸気時の胸部ＣＴ画像（第２ボリュームデータ）と呼気時の胸部ＣＴ画像（第１ボリュームデータ）とが記憶装置１１２を介して画像位置合わせ処理部１１９に送られ、吸気時の胸部ＣＴ画像（第１ボリュームデータ）に呼気時の胸部ＣＴ画像（第２ボリュームデータ）を位置合わせする（Ｓ３）。この位置合わせには両画像間で呼吸に伴う肺野変形が存在するので、非剛体位置合わせ（non-rigid registration）の技術を用いてボクセルレベルで位置合わせを行う。この位置合わせ技術は肺の領域抽出、気管支の領域抽出、あるいは肺動脈の領域抽出を行い、それらの領域の呼気時、吸気時画像の差から呼気時画像空間のワープフィールドを作成して呼気時画像を変形する手法などを用いる。

また吸気時の胸部ＣＴ画像（第２ボリュームデータ）が記憶装置１１２を介して領域分割処理部１１７に送られ、肺野領域が抽出され、肺野領域が左肺領域と右肺領域とに分割され、左右肺領域それぞれはさらに図５に例示するように右肺は上葉領域、中葉領域、下葉領域の３つの肺葉領域、左肺は上葉領域、下葉領域の２つの肺葉領域に分割され、さらに葉気管支の分岐点により肺葉領域はさらに小さな単位である区域領域に分割される。

吸気時の胸部ＣＴ画像（第２ボリュームデータ）と呼気時の胸部ＣＴ画像（第１ボリュームデータ）とが記憶装置１１２を介して画像位置合わせ処理部１１９に送られ、吸気時の胸部ＣＴ画像（第２ボリュームデータ）に呼気時の胸部ＣＴ画像（第１ボリュームデータ）を位置合わせする（Ｓ３）。この位置合わせには両画像間で呼吸に伴う肺野変形が存在するので、非剛体位置合わせ（non-rigid registration）の技術を用いてボクセルレベルで位置合わせを行う。この位置合わせ技術は肺の領域抽出、気管支の領域抽出、あるいは肺動脈の領域抽出を行い、それらの領域の呼気時、吸気時画像の差から呼気時画像空間のワープフィールドを作成して呼気時画像を変形する手法などを用いる。

次に、後述する工程Ｖ３，Ｖ４において利用されるボクセルの動き量（動きベクトルの距離）に関連する分析ではそれぞれの肺構造の位置において適正な重症度を示すテーブル（以下、重症度テーブルと記す）を必要とするため、ボクセル単位、あるいは例えば４×４×４ボクセルによりなるブロック領域単位でその適正値を有する形状・構造の標準胸部ＣＴ画像を吸気時胸部ＣＴ画像に位置合わせを行う。位置合わせ手法は工程Ｓ３と同様の処理が適用される。

以降は処理内容をボクセル単位のカラー画像を生成する工程Ｖ１−Ｖ４と領域単位のカラー画像を生成する工程Ｒ１−Ｒ４に分けて個々に説明する。各々の処理概要は以下の通りである。即ち、ボクセル単位では各ボクセルの動きベクトル量（変位距離）と、ＣＴ値の分析から得た一つの特徴値と、体積変化率との少なくとも２つの組み合わせから、工程Ｓ５で得られた標準画像のボクセルに対応した予め定められている重症度テーブル上にマッピングして、その画素に重症度に応じたカラーを割り当ててカラー表示する。一方、領域単位では領域としての動きベクトルと個々の領域で計算できるＣＴ値を利用した低吸収領域率などの評価値の分析から得た一つの値を工程Ｓ４で得られた分割領域毎又は工程Ｓ５で得られた標準画像のブロック領域に対応した重症度テーブル上にマッピングして、その領域に重症度に応じたカラーを割り当ててカラー表示する。

上記において予め定められている重症度テーブルは、それぞれ予め得られた多数のGOLD評価がなされている臨床データに対して同一の処理（S1-5, V1-2, R1-2）を行い、それぞれ動きやCT値などのパラメータの組み合わせに応じたGOLD評価との主成分分析を行なって作成される。

まずボクセル単位を処理対象とした工程Ｖ１では、位置合わせされた２つの胸部ＣＴ画像からベクトル計算処理部１２１によりボクセル（画素）レベルの動きベクトルを計算し、その３つの方向成分であるＸ，Ｙ，Ｚ方向の動き量及びベクトル絶対値を計算する。すなわち、動きベクトルは、例えば、横隔膜移動方向を示すＺ軸方向成分と、胸骨移動方向を示すＹ軸方向成分と、前記Ｚ軸及びＹ軸に直交するＸ軸方向成分とに分化される。

次に、工程Ｖ２では、位置合わせされた２つの胸部ＣＴ画像から重症度レベル計算部１２３でボクセル（画素）単位の吸気時胸部ＣＴ画像及び呼気時胸部ＣＴ画像のＣＴ値とその差分値を計算する。また第１、第２のボリュームデータ間での領域の体積変化率を計算する。

工程Ｖ３では重症度レベル決定部１２３によりボクセル単位で重症度レベルが決定され、それに応じたカラー画像が作成され表示される。この工程Ｖ３では、ボクセル毎の動きベクトルで示された変位方向に係る変位距離と、動きベクトルに関するＸＹＺ直交３軸に関する３つの方向成分（距離）との中の２つの方向成分に基づいて、呼吸に伴う各部の移動との観点からボクセル毎に重症度レベルが前記の予め得られている重症度テーブル上にそのボクセルの分析値をマッピングすることで決定（分類）される。例えば、図６に示すように、Ｙ方向の動き量（移動距離）とＺ方向の動き量とからＧｏｌｄ０乃至Ｇｏｌｄ３の４つの重症度レベル（ステージ）がボクセル毎に決定される。例えばＧｏｌｄ０からＧｏｌｄ３に向かって重症度が高くなる。

図７に示すように、Ｘ方向の動き量とＺ方向の動き量とから重症度レベルを決定してもよいし、図８に示すように、ＸＹＺの一方向の動き量と動きベクトル方向に係る動き量とから重症度レベルを決定してもよい。

ボクセル毎に決定された重症度レベルを第一の吸気時胸部CT画像上に重症度レベルに応じたカラーを割り当ててカラー画像を作成し、図１４に示すように表示する。この重症度カラー割り当ては重症度レベルＧｏｌｄ３が例えば赤色、Ｇｏｌｄ０には緑色等４つの重症度レベルに異なるカラーを割り当てる。またはＧｏｌｄ３からＧｏｌｄ０に向かって連続的なグラデーションとして例えば赤から緑に変化させるものであっても良い。尚、この例ではGold 0〜3の４段階としているが、Gold 0〜4の５段階で領域分けを行なっても良い。

上述の説明では、ベクトルで示された変位距離と動きベクトルに関するＸＹＺ直交３軸に関する３つの方向成分の中から２つのパラメータに基づいて、閉塞性肺疾患の重症度に関する重症度レベルを決定したが、動きベクトル量等の中から３つのパラメータに基づいて重症度レベルを決定してもよい。

工程Ｖ４では、各ボクセルの動きにおいて上記と同様に動きベクトル方向に関する動きベクトル量、動きベクトルに関するＸＹＺ直交３軸に関する３つの方向成分の中から選択された１つ以上の動き方向成分と、ＣＴ値から求められた特徴値とを用いて重症度レベルが決定される。例えば、図９に示すように、ＸＹＺの一方向の動き量若しくは動きベクトル方向に係る動き量と、２つの胸部ＣＴ画像の一方のＣＴ値若しくはその差値とから重症度レベルを決定してもよい。図１０ではz方向の動き量に対するｙ方向の動き量にＣＴ値又はその２画像間でのＣＴ値差を３軸目としてそれら３つのパラメータから決まる４領域に従ってレベルを決定する。図１１ではz方向の動き量に対するｘ方向の動き量にＣＴ値を３軸目としてもよいし、図１２ではボクセルの動きベクトル方向の動き量（即ち、動きベクトルの絶対値）に対するｘ，ｙ，ｚ方向の一つの方向成分のなかの最大量にＣＴ値を３軸目としてもよい。

こうして２つ又は３つのパラメータからボクセル毎に決定された重症度レベルを第一の吸気時胸部CT画像上に重症度レベルに応じたカラーを割り当ててカラー画像を作成し、図１４に示すように表示する。このカラー設定はGold 3が例えば赤色からGold 0の緑色と４つのGold領域に単独にカラーを定めるか、あるいはGold 3の３軸交点からGold 0の手前に連続的に赤から緑に変えても良い。

上記の３つのパラメータに加えて、呼気時のＣＴ値などを加えて、４つのパラメータ、さらに他のパラメータを加えて５次元以上の多次元空間の中で主成分分析などの手法を用いてGold値を設定してもよい。そのGold値に対応したカラーを当該ボクセルのカラー値とする。

次に領域単位でレベルを決定してカラー画像を生成する処理について説明する。工程Ｓ４で分割した分割領域、あるいはブロック領域に関する構造体体積（またはその変化）、低線量領域の割合を表すLAV%値（またはその変化）などのパラメータを、重症度レベルの決定のための選択肢に加える。

工程Ｒ１では、Ｓ４で得られた左右肺、肺葉又は区域の分割領域又はブロック領域に関する各領域としての動きベクトルを計算する。この計算は上記の各領域内の複数のボクセルに関する複数のベクトルの動きベクトルの平均値、中間値などを当該領域に関する動きベクトルとして捉える。工程Ｒ２では、Ｓ４で得られた左右肺、肺葉又は区域に関する各領域毎又はブロック領域毎のＣＴ値から各領域における低線量領域の割合を表す吸気時LAV%値（閾値を例えば-950HU）、呼気時LAV%（閾値を例えば-856HU）を求める。あるいは各領域における吸気時、呼気時のCT値の平均値、分散、尖度などを評価値として領域評価値計算部１２０により計算する。

工程Ｒ３は、工程Ｖ３と同様に、重症度レベル決定部１２３により領域単位で重症度レベルが決定され、それに応じたカラー画像が作成され、図１５に例示するように表示される。この工程Ｒ３では、領域毎の動きベクトルの絶対値と、領域毎の動きベクトルに関するＸＹＺ直交３軸に関する３つの方向成分（距離）との複数の動きパラメータの中から２つの方向成分に基づいて、領域毎に重症度レベルが決定される。

工程Ｒ４では、Ｒ３の複数の動きパラメータの中から１つ以上の方向成分に、例えば領域毎の吸気時ＣＴ値平均値、その差、LAV%値から選択したパラメータを加えて、これら２つ以上のパラメータから重症度レベルが決定される。このように多次元空間上で決定した重症度レベルに応じてカラー画像を作成し、表示するようにしてもよい。

このR3,R4の領域単位での２次元カラー画像、３次元カラー画像は上記ボクセルまたはブロック単位の重症度テーブルとは異なり、それぞれの肺葉や区域などの分割領域に対応する重症度テーブルを用いることになる。

動き、ＣＴ値、構造物など複数のパラメータを用いることにより、より正確な局所での慢性閉塞性肺疾患の診断評価が可能となる。

上述の説明では、ボクセル単位、ブロック領域単位、分割領域単位で重症度レベルが決定され、その空間分布をカラー画像として生成したが、肺全体として重症度レベルを決定するようにしてもよい。肺全体としての重症度レベルは、領域単位の重症度レベルの決定方法と同じ処理により決定され得る。特に肺全体の場合はカラー表示だけでなく単にGOLD値を提示することや本手法による肺全体の分析値を点として２次元又は３次元の重症度テーブルの場合はそのテーブル上に表示しても良い。

図１８に示すように、領域の体積変化率と、ＣＴ値から求めた特徴値とから重症度レベルを領域ごとに決定してもよい。図１９に示すように、領域の体積変化率と、領域ごとの動きベクトルの距離（動き量）と、領域ごとの動きベクトルのＸＹＺの中の一方向成分とから重症度レベルを領域ごとに決定してもよい。

なお、上述の説明では、呼吸時相の異なる２つの画像間で動き量等を求めて、重症度レベルを決定していたが、図２０に示すように、呼吸時相Ｖt1、Ｖt2、Ｖt3の異なる３つの画像、さらに多くの画像間で動き量等を求めて、重症度レベルを決定してもよい。呼吸時相の異なる３つの画像から領域ごとの動きベクトルに関するＸＹＺの一方向成分の動き量は、ＴＲ１、ＴＲ２，ＴＲ３のいずれかが選択される。例えば動きベクトルに関するＸ方向成分の動き量は、呼吸時相Ｖt1の画像上のあるボクセルのＸ座標をＸ１，呼吸時相Ｖt2の画像上の対応ボクセルのＸ座標をＸ２，呼吸時相Ｖt3の画像上の対応ボクセルのＸ座標をＸ３とすると、
ＴＲ１＝｜Ｘ１−Ｘ２｜＋｜Ｘ２−Ｘ３｜
ＴＲ２＝｜Ｘ１−Ｘ２｜／｜Ｘ２−Ｘ３｜
ＴＲ３＝（｜Ｘ１−Ｘ２｜＋｜Ｘ２−Ｘ３｜）／２
で与えられる。つまり、複数の画像間の移動距離の総和、複数の画像間の移動距離の比率、複数の画像間の移動距離の平均で重症度レベルを決定するパラメータとしての動き量が求められる。

本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれると同様に、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。

１００…架台部（ガントリ）、１０１…Ｘ線管、１０２…回転リング、１０３…Ｘ線検出器、１０４…データ収集回路、１０５…非接触データ伝送装置、１０６…前処理部、１０７…架台駆動部、１０８…スリップリング、１０９…高電圧発生部、１１０…ホストコントローラ、１１２…記憶装置、１１５…入力装置、１１９…画像位置合わせ処理部、１１７…領域分割処理部、１２０…領域評価値計算部、１２１…ベクトル計算処理部、１２２…表示部、１２３…重症度レベル決定部、１２４…カラー画像生成処理部。

Claims (25)

1. 被検体胸部の形態を表す呼吸時相の異なる複数の画像のデータを記憶する画像記憶部と、
   前記複数の画像間での部位の動き量を画素又は領域ごとに計算する動き量計算部と、
   前記動き量と、前記画像の画素値から求められる特徴値と、前記領域の大きさの変化率との少なくとも２つの組み合わせにより慢性閉塞性肺疾患の重症度に関するレベルを画素又は領域ごとに決定するレベル決定部と、
   前記決定されたレベルに関する情報を出力する出力部とを具備することを特徴とする医用画像処理装置。
2. 前記特徴値は、前記画素値、前記画素値の差、前記領域内の画素数に対する閾値未満の画素値を有する画素の数の割合、前記領域内の画素値分布の分散、または前記領域内の画素値の平均値であることを特徴とする請求項１記載の医用画像処理装置。
3. 前記決定されたレベルに関する２次元又は３次元のマップを発生するマップ発生部をさらに備えることを特徴とする請求項１記載の医用画像処理装置。
4. 前記画像は２次元画像又は３次元画像であることを特徴とする請求項１記載の医用画像処理装置。
5. 前記動き量計算部は、前記画素又は領域ごとにベクトルを計算することを特徴とする請求項１記載の医用画像処理装置。
6. 前記動き量は、前記ベクトルの一方向成分又は二方向成分で表されることを特徴とする請求項５記載の医用画像処理装置。
7. 前記領域の大きさの変化率は、前記複数の画像内の一の画像における近傍する複数の画素で区画される領域の面積又は体積に対する、前記近傍する複数の画素それぞれに前記ベクトルにより対応する他の画像における複数の画素で区画される領域の面積又は体積の比率により計算されることを特徴とする請求項５記載の医用画像処理装置。
8. 前記複数の画像には吸気時に発生された吸気画像と呼気時に発生された呼気画像とが含まれ、
   前記動き量、前記特徴値及び前記変化率は前記吸気画像と前記呼気画像とを用いて計算されることを特徴とする請求項１記載の医用画像処理装置。
9. 前記領域として肺野の右肺領域と左肺領域とが適用されることを特徴とする請求項１記載の医用画像処理装置。
10. 前記右肺領域と前記左肺領域それぞれは複数の肺葉領域に分割されることを特徴とする請求項９記載の医用画像処理装置。
11. 前記肺葉領域それぞれは複数の肺区域領域に分割されることを特徴とする請求項１０記載の医用画像処理装置。
12. 前記領域として近傍する所定数のボクセルからなる集合体としてのブロック領域が適用されることを特徴とする請求項１記載の医用画像処理装置。
13. 前記複数の画像間で解剖学上の位置を合わせる位置合わせ処理部をさらに備え、
    前記動き量計算部は、前記位置合わせされた画像間での部位の動き量を計算することを特徴とする請求項１記載の医用画像処理装置。
14. 被検体胸部の形態を表す呼吸時相の異なる第１，第２のボリュームデータを記憶する記憶部と、
    前記第１，第２のボリュームデータ間での部位の動きベクトルを画素又は領域ごとに計算するベクトル計算部と、
    前記動きベクトルと、前記第１，第２のボリュームデータの画素値から求められる特徴値と、前記領域の大きさの変化率との少なくとも２つの組み合わせにより慢性閉塞性肺疾患の重症度に関するレベルを画素又は領域ごとに決定するレベル決定部と、
    前記決定されたレベルに関する空間的な分布を示すマップを発生するマップ発生部とを具備することを特徴とする医用画像処理装置。
15. 前記第１のボリュームデータは吸気時に収集されたものであり、前記第２のボリュームデータは呼気時に収集されたものであることを特徴とする請求項１４記載の医用画像処理装置。
16. 前記レベル決定部は、前記特徴値と、前記動きベクトルに関する一方向成分とに基づいて前記レベルを決定することを特徴とする請求項１４記載の医用画像処理装置。
17. 前記レベル決定部は、前記特徴値と、前記動きベクトルに関する二方向成分とに基づいて前記レベルを決定することを特徴とする請求項１４記載の医用画像処理装置。
18. 前記特徴値は、前記画素値、前記画素値の差、前記領域内の画素数に対する閾値未満の画素値を有する画素の数の割合、前記領域内の画素値分布の分散、または前記領域内の画素値の平均値であることを特徴とする請求項１４記載の医用画像処理装置。
19. 前記動きベクトルは、前記領域各々に含まれる複数のベクトルから、平均距離を示すベクトル、距離の中間値を示すベクトル又は距離の中央値を示すベクトルが選択されることを特徴とする請求項１４記載の医用画像処理装置。
20. 前記第１、第２のボリュームデータは、Ｘ線コンピュータ断層撮影装置又は磁気共鳴イメージング装置により発生されることを特徴とする請求項１４記載の医用画像処理装置。
21. 前記領域として肺野領域は右肺領域と左肺領域とに分割され、前記右肺領域と前記左肺領域それぞれは複数の肺葉領域に分割され、または前記肺葉領域それぞれは複数の肺区域領域に分割されることを特徴とする請求項１４記載の医用画像処理装置。
22. 前記動きベクトルは、横隔膜移動方向を示すＺ軸方向成分と、胸骨移動方向を示すＹ軸方向成分と、前記Ｚ軸及びＹ軸に直交するＸ軸方向成分とに分化されることを特徴とする請求項１４記載の医用画像処理装置。
23. 前記第１，第２のボリュームデータ間で解剖学上の位置を合わせる位置合わせ処理部をさらに備え、
    前記ベクトル計算部は、前記位置合わせされた画像間での部位の動きベクトルを計算することを特徴とする請求項１４記載の医用画像処理装置。
24. 前記領域の大きさの変化率は、前記第１のボリュームデータにおける近傍する複数のボクセルで区画される領域の面積又は体積に対する、前記近傍する複数のボクセルそれぞれに前記動きベクトルにより対応する前記第２のボリュームデータにおける複数のボクセルで区画される領域の面積又は体積の比率により計算されることを特徴とする請求項１４記載の医用画像処理装置。
25. 被検体胸部の形態を表す呼吸時相の異なる第１、第２のボリュームデータ間で解剖学上の位置を合わせ、
    前記位置を合わされた第１、第２のボリュームデータ間でボクセル単位又は領域単位で動きベクトルを計算し、
    前記動きベクトルと、前記第１，第２のボリュームデータの画素値から求められる特徴値と、前記領域の大きさの変化率との少なくとも２つの組み合わせにより慢性閉塞性肺疾患の重症度に関するレベルを画素又は領域ごとに決定し、
    前記決定されたレベルに関する空間的な分布を示すマップを発生する
    ことを特徴とする医用画像処理方法。