JP4572294B2 - 画像処理プログラム及び画像処理方法 - Google Patents

Description

本発明は、３次元ＣＴ装置や３次元ＭＲＩ装置等の断層撮像装置を用いて得られた画像データの処理プログラム及び画像の処理方法に関する。

人体の断層像を取得する断層撮像装置としてはＸ線ＣＴ装置、ＭＲＩ装置が広く普及しており、特に近年では人体の断層像を複数重ね合わせ３次元画像データを得ることのできる３次元Ｘ線ＣＴ装置、３次元ＭＲＩ装置が実現している（例えば下記非特許文献１、２参照）。しかし、従来の３次元ＣＴ装置、３次元ＭＲＩ装置では撮像対象領域が狭く、しかもその撮影に数秒以上の時間を要しており、動きのある被写体の撮影は困難であった。

これに対し最近になって、広い範囲の撮像が可能でかつ高速連続撮影可能な３次元ＣＴ装置の研究開発が行われ、試作機が登場してきた。この装置によって、人体を動かしながら撮影を行い、対象領域の３次元画像データを刻々と取得できるようになった。なお、本明細書では時系列的に得られる３次元画像データを４次元画像データといい、この画像を得ることができる断層撮像装置を４次元断層撮像装置という。特にＣＴ装置の場合は４次元ＣＴ装置、ＭＲＩ装置の場合は４次元ＭＲＩ装置のように呼ぶこととする。
飯沼武，舘野 之男編著 「Ｘ線イメージング」コロナ社 ２００１年 レイ・Ｈ．ハシェミ他著 「ＭＲＩの基本パワーテキスト」メディカルサイエンスインターナショナル社 １９９９年

これら装置により得られた画像から被写体の構成要素の動態解析を行うことで、従来困難であった医学的知見が得られる可能性がある。このためには、各構成要素を抽出する技術およびそれら時間変化する構成要素を対応付ける技術が必要であるが、新規に登場した４次元断層撮像装置に対してそれらの技術開発は全くなされていない。

そこで、本発明は上記課題を解決し、複数の３次元画像データの処理を行い被写体の動態解析を行うことが可能となる画像処理プログラム、画像処理方法を提供することを目的とする。

上記目的を達成するために、本発明は具体的には以下の手段を採用する。
まず、コンピュータに、時系列的に得られた複数の３次元画像データのそれぞれにおいて抽出対象領域を分割するステップ、分割された複数の前記抽出対象領域の変位量を求めるステップと、を実行させるための画像処理プログラム、とする。

なお、本手段において、時系列的に得られた複数の３次元画像データのそれぞれにおいて抽出対象領域を分割するステップは、３次元画像データにおける最大の画素値を求めるステップ、求めた最大の画素値よりも小さい所定の画素値を設定するステップ、３次元画像データを複数の２次元画像データに分割してそのいずれかを選択し、選択された２次元画像データのいずれか一側辺から順次走査して所定の画素値以上となった画素を開始点と定め、開始点となる画素と隣接する画素における画素値が上述の所定の画素値以上であると判断した場合、更に画素に隣接する画素の画素値に対し同様に判断を繰返すステップ、を有することも望ましい。

また、本手段において、変位量は、回転及び平行移動の量に基づいて求められること、変位量は、非線形変形のパラメータにも基づいて求められることも望ましい。

また、本手段において、分割された前記抽出対象領域の変位量を求めるステップは、時系列的に得られた複数の３次元画像データそれぞれに対し、分割された抽出対象領域をスライスする２次元画像データを複数選択し、そのそれぞれの２次元画像データにおける抽出対象領域の重心点を求め、求めた複数の重心点から抽出対象領域の軸を定め、軸の変位量を求めるステップ、を有することも望ましい。

また抽出対象領域は３次元画像データのそれぞれにおいて複数抽出されることも望ましく、その際に、分割された複数の抽出対象領域の変位量を求めるステップは、時系列的に得られた複数の３次元画像データそれぞれに対し、分割された抽出対象領域をスライスする２次元画像データを複数選択し、そのそれぞれの２次元画像データにおける抽出対象領域の重心点を求め、求めた複数の重心点から抽出対象領域の軸を定め、軸の変位量を求めるステップ、を有すること、分割された複数の前記抽出対象領域の時系列的な変位量を求めるステップは、更に、時系列的に得られた複数の３次元画像データにおいて分割された抽出対象領域同士の最短距離を求めるステップ、を有することも望ましい。

以上により本発明は、複数の３次元画像データの処理を行い被写体の動態解析を行うことが可能となる画像処理プログラム、画像処理方法を提供することができる。

以下、本発明の実施の形態の一例について、図面を用いて説明する。
本実施形態に係る画像処理プログラムは、時系列的に得られる３次元画像データ（４次元画像データ）に対し、各時刻における３次元画像データ内の構成要素を抽出し、それらの形状の時間変化を追うことで構成要素の形状変化を取得できるものである。この方法の適用先としては、特に最近研究開発がなされ試作機が作られている４次元ＣＴ装置が好適であり、本実施形態では４次元ＣＴ装置に適用した例を具体的に説明する。

まず４次元ＣＴ装置の基本原理と試作機での性能を説明する。
図1は４次元ＣＴ装置の原理を示すための図である。図１（ａ）で示すように、本実施形態に係る４次元ＣＴ装置は、面状にＸ線を放射するＸ線源１と、このＸ線源が放射するＸ線を検出する２次元の検出器２と、を有して構成されており、Ｘ線源１と２次元検出器２とは測定対象となる被写体３を間に挟んで回転する（このように２次元検出器とＸ線管を被写体の周りに回転させる方式を採用するＣＴ装置は一般にコーンビームＣＴ装置とも呼ばれている）。より具体的に説明すると、図１（ｂ）で示すように、本実施形態に係る４次元ＣＴ装置はＸ線源１と２次元検出装置２とがリング４に固定されており、リング４を回転させてもそれらの位置関係を一定に保つことができる。そして一回の回転中複数の２次元のＸ線強度データを取得し、それらに対し画像再構成処理を行うことで３次元ＣＴ画像データを得ることができる。被写体３は測定中においてリング４の回転中心近傍で運動を行うが、運動は時間分解されそれぞれの時間における３次元画像データとして複数取得され、４次元画像データとなる。なお、２次元のＸ線強度データ、３次元画像データ、４次元画像データなどの各種画像データ、更にはこれら各種画像データを処理するためのプログラムは断層撮像装置内の記録媒体若しくは断層撮像装置に電気的に接続されるパーソナルコンピュータの記録媒体に格納され、必要に応じて加工又は実行され、各種画像データ、プログラム以外の部分については周知の構成を採用することができる。

なおここで、図２に各３次元ＣＴ画像データ構成について図を用いて説明する。図２で示すように、本実施形態に係る３次元ＣＴ画像データは複数の２次元ＣＴ画像データの積み重ねとして表現されたものとなっている。なお２次元画像データとは２次元的に配列した複数の画素からなる画像データであり、３次元画像データとは３次元的に配列した複数の画素からなる画像データである。なお図２は人間の膝関節における３次元ＣＴ画像データ構成の模式図である。

本実施例に係る４次元ＣＴ装置では、ドラムの回転を高速にしかも複数回回転させ、それぞれの回転における３次元画像を取得し、時系列に並んだ３次元ＣＴ画像データ即ち４次元ＣＴ画像データを得ることができるため、これらが更に所定の時間に相当する分だけ取得されることとなる。本実施形態に係る４次元ＣＴ装置では、幅２３ｃｍ、（ドラム円周方向の）長さ９０ｃｍの２次元検出器を回転数２回／秒回転させ、サイズが回転面内で約２４０ｍｍ×２４０ｍｍ（５１２×５１２画素）、回転軸方向に１２８ｍｍ（２５６スライス）、空間分解能等方的に約０．５ｍｍ、時間分解能０．５秒の４次元ＣＴ画像データを得ることができている。なおドラムの直径や検出器のサイズ等といった諸条件については４次元ＣＴ画像データをとることができる限りにおいて特段制限がされるわけではなく、適宜調整が可能である。

以上、被写体が動いている場合（例えば膝関節が運動している場合）であっても撮影することが可能となり、その時系列に並んだ３次元ＣＴ画像データ（４次元ＣＴ画像データ）データが得られる。

次に、この４次元ＣＴ画像データを処理する本実施形態に係る画像処理プログラム（以下「本プログラム」ともいう。）の実行による処理について具体的に説明する。

本プログラムの処理は、大きく２つのステップを有して構成される。第１のステップは各フレーム（各時点）における３次元画像データから分析対象となる構成要素を抽出するステップであり、第２のステップは抽出結果を用いて種々の解析を行うステップである。なお図３に本プログラムのステップに関するフローチャートの概要を示し、詳細について以下説明する。本プログラムはこのように構成要素を抽出し、骨や血液等、着目したい構成要件のみについて時系列データを取得し、詳細に検討することを可能とする。

まず、各構成要素を抽出する第1のステップについて述べる。本ステップは、対象とする構成要素が有する画素値の範囲と、空間位置の特徴に注目し、この要素を簡便に抽出できるよう工夫している。ここで対象とする構成要素としては例えば膝関節等の骨、血管、神経等人体の様々な部位が該当し、これらは３次元画像データとして把握されたとき、３次元的に連結された領域として表される。なお、ここで画素値とは、Ｘ線の吸収係数に対応した値をいい、Ｘ線の吸収係数が大きいほど大きな値として把握されるものである。

上記３次元的に連結された領域の抽出法としては、一般的な領域拡張法を基礎として利用することができる。まず、一般的な領域拡張法について図４を用いて説明する。

この方法では、まず適当な開始点を抽出したい構成要素の領域（以下「抽出対象領域」という。）の一点として選択し、その次に開始点と隣接している画素に着目し、その隣接する画素の値が所定の画素値の範囲にある場合は抽出対象領域として、範囲にない場合は非抽出対象領域とする。そして対象領域として抽出された場合は更にその隣接する画素の値について先ほどと同様に所定の画素値の範囲か否かを判断することを繰返し、対象領域を拡張させていく。これにより、抽出対象領域を特定する。

より具体的に説明する。図４で示すように、ある位置（ｘ、ｙ、ｚ）を開始点５とし、その点における画素値ｆ（ｘ、ｙ、ｚ）が下記式（１）を満たす場合、当該位置が抽出対象領域内の一点であると判断する。

そして、上記式（１）を満たしていると判断された場合、Ｌ（ｘ、ｙ、ｚ）のラベルを付す。Ｌ（ｘ、ｙ、ｚ）は、この点が抽出対象領域であることを意味するラベルである。そして、このラベルが付された座標に隣接する画素それぞれに対し、下記式（２）で示される判断処理を行いこれを繰り返すことにより、抽出対象領域を定めていくことができるのである。

この方法では回始点の設定を手動で行うのが通常である。しかしながら、本実施形態に係る４次元ＣＴ画像データは、３次元ＣＴ画像データが更に時系列で多数並べられたものであるため、時々刻々対象領域が変化する各３次元ＣＴ画像データ各々に対し手動で開始点を定めていくことは極めて困難である。そこで本実施形態では、以下の方法により自動化を可能とした。図５に本実施形態にかかる開始点探索のアルゴリズムを説明するための模式図を示す。

本プログラムの第１のステップでは、まず４次元ＣＴ画像データのうちのいずれかの３次元ＣＴ画像をひとつ選択し、その３次元画像データのうち最も大きな画素値を探す。そしてこの画素値を求め、この画素値に１以下の適切な範囲の値を設定し、抽出対象領域における画素値の範囲を定める。これにより上記式（１）及び（２）の判断の範囲を定めることができる。

次に、上記画素値の範囲を定めた後、開始点の探索を開始する。本ステップでは３次元ＣＴ画像データを構成する複数の２次元ＣＴ画像データひとつを選択し、更にその２ＣＴ次元画像データの一辺の画素列又は行から順次開始点の捜査を行う。具体的に図５のイメージ図で説明すると、まず３ＣＴ次元画像データを構成する複数の２次元画像データのうちいずれかの一側面に相当する２次ＣＴ元画像データ６を選択する。そしてその２次ＣＴ元画像データ６のうち、適当な一辺の行７を選択し、その行７において所望の画素値の条件を満たす開始点を順次探索し（図中（１）の方向）、条件を満たす点が無い場合は次の行へと移動し（図中（２）の方向）、開始点を繰り返し探索していく。なお、それでも２次元ＣＴ画像データ６において予め定めた開始点の条件を満たした点が見つからない場合は、隣の２次元ＣＴ画像データに移り（図中（３）の方向）、開始点を探索し、見つかるまでこれを繰り返す。これにより、高い精度で抽出対象領域を定めることができる。

３次元ＣＴ画像データを構成する２次元ＣＴ画像データのいずれを選択するか、２次元ＣＴ画像データのうちのどの一辺の画素列又は行を選択するかについては、抽出対象領域としたい構成要素が３次元ＣＴ画像データのどの位置にあるか３次元画像データを視覚化し予め定めておくことが有用である。なお図５の例は、大腿骨と脛骨とが描出されている膝関節近傍の３次元ＣＴ画像データであって（但し大腿骨については見えない）、符号６で示される２次元ＣＴ画像は脛骨を下側から見た場合の２次元ＣＴ画像データである。

即ち本プログラムはこのような処理を行うことで開始点が容易に探索可能となり、更に３次元の領域拡張法を適用することで非常に煩雑な手動による開始点の指定の繰り返しを必要とすることなく簡便に脛骨の抽出を行うことができる。なお、大体骨や膝蓋骨など他の骨についても同様の考え方に基づき、適切な２次元ＣＴ画像データ、およびその一辺を選択し、そこからスキャンを開始し、高い精度で大体骨における開始点を見つけることができる。

なお上述のように本実施形態では予め３次元画像データを視覚化し、開始点探索を行う２次元ＣＴ画像データの一辺を定めておくことが有用であるが、一つの抽出対象領域に該当する骨に対して探索を開始する２次元ＣＴ画像データの一辺を選択させ、一つの抽出対象領域を決定した後、この抽出対象領域を抽出対象から外し、更に同じ一辺から開始点の探索を行うようにすることが可能であり、このようにすれば一度選択した抽出対象領域を再度抽出対象領域として選んでしまうことを防ぎ、これ以外の他の骨を抽出対象領域として選定することが非常に容易にできるようになる。特に所望の範囲に含まれる画素値を示す構成要素を抽出対象領域としたい場合であれば、２次元ＣＴ画像データの一側辺を選択する作業さえも省略可能となる。なお図６に膝関節における抽出対象領域を分割した結果について示す。図６は大腿骨および脛骨についての分割を環状断（ｃｏｒｏｎａｌ）面から見た場合、矢状断（ｓａｇｉｔｔａｌ）面から見た場合の両方を示している。また図７にこの領域分割の結果から得られる表面画像を示す。なお、本プログラムは各時刻における表画像を時系列的に表示させていくことにより、動画像として表示させることができる。

なお本ステップは高いＸ線の吸収を示す骨について行うことで画像データの取得に特段の処理を行うことなく開始点を定めることができたが、Ｘ線の吸収のあまり高くない部位に対しても同様に行いたい場合、例えば血液や関節軟骨等に対しても行いたい場合は、造影剤等の投与を行うことにより血液や関節軟骨の画素値を高めることができ、本ステップを適用させることができる。

次に第２のステップについて説明する。本ステップでは多様な解析が可能であるが，まず基本になるのがフレーム間で対応する領域ごとに位置合わせを行い、分割された抽出対象領域の時間的な変位量（以下単に「変位量」ともいう。）を算出することである。まずここでは運動の前後において局所的な変形がないとみなせる骨を対象とし、剛体としての回転及び平行移動で変位量を算出する。

ここで図８を用いて変位量の算出について説明する。図８は変位量の算出の原理を示す概念図である。まず移動量算出にあたり、移動前の画像と移動後の画像とを、それぞれ、f₁(p₁)、f₂(p₂)と表す。ここでp_1、p₂は画素位置を表す列ベクトルであり、ｘｙｚを用いて表すと、p₁=[x₁、y₁、z₁]^ｔp2=[ x₂、y₂、z₂]^t、と表すことができる。ただし、[ ]^tは行列、ベクトルの転置を表す。

上記の表現の下、剛体変換は次式で与えられる。

なおＲは回転を表す３×３の行列であって、ｚ軸、ｘ軸、ｙ軸のまわりの回転をそれぞれθ、α、βとして次のような式で表される。

また，tは移動を表すベクトルであり，ｘ軸、ｙ軸、ｚ軸各方向の移動量を用いて次のような式で表される。

以上、上記の剛体変形によって移動前の画像が変換され、移動後の画像に一致したとすれば、その場合それら画素毎の画素値の差の二乗和、即ち
が０となるはずである。しかし実際には、上記の評価関数を最小とするように、回転と平行移動のパラメータθ、α、β、ｔ_ｘ、ｔ_ｙ、ｔ_ｚの６つを適当な最適化のアルゴリズムを用いて決定する。なお、上記の評価関数を最小とするかわりに、２枚の画像の相互相関を最大にするようにパラメータを求めてもよい。以上の処理により、２つの物体間の変位量が得られる。

以上の基本処理が終わった後には、医学的に有用な種々の解析が可能になる。

本実施形態では、抽出対象領域として剛体と考えることができる骨を用いている。即ち骨の場合は回転と平行移動のみで変位量を容易に求めることができる。しかし一方、血管や軟骨などの軟体の場合は局所的な変形が入ってくるため、回転と平行移動だけでは正確な変位量を求めることは難しいが、上記に加え、非線形変換のパラメータを更に加えることで変位量を得ることができる。具体的には、対象領域上の各点に対して変位量を求めることで実現できる。以下に説明する。

まずこの場合、位置p₁における変形量のｘ、ｙ、ｚ成分を
と書くことにすると、剛体変換の際の変換式p₂=Rp₁+tに対応する式は
と書ける。また，剛体変換の際の評価関数は
と書き直すことができる。ここで、求めるべきパラメータはΔp(p₁)となる（図９参照）。これにより剛体だけでなく軟体に対しても変位量を求めることができるようになる。

なお、この方法で局所的な変形の情報が得られた場合、予め軟骨等がもつ弾性係数等の力学的諸係数を別途計測し、これら諸係数の分布と構成要素の各点に加わる力の分布を与えることで変形量の分布を微分方程式にしたがって決定することができ、逆に、上記方法で得られた変形量の分布から微分方程式を解くことによって構成要素の各点に加わっている力のベクトル分布を推定することが可能となる（図１０参照）。

次に、膝の屈曲角に対する大腿骨−脛骨間の最近接点間距離および最近接点位置の推移をグラフ化する方法について図１１を用いて説明する。

膝の屈曲角を得るためには、まず、大腿骨と脛骨の骨軸を取得する必要がある。骨軸抽出は以下のような方法で行う。まず、抽出対象領域である骨領域をスライスする任意の枚数（Ｎ枚）の２次元ＣＴ画像データより、骨領域の重心点を数学的に算出し、求めた複数の重心点の３次元的なばらつきに対して主成分分析を行い、結果として得られた第一主成分の軸を骨軸とする（図１１（ａ））。この軸は大腿骨と脛骨、関係を調べたい複数の骨に対して行う。そして上記の方法を骨軸に対して同様に用いて得られた変位量を用いてこの骨軸を剛体変換し例えば膝の屈曲角θを得ることができる。これにより屈曲角による膝の解析が可能となる（図１１（ｂ））。

また次に、大腿骨-脛骨間の最短距離を外側顆、内側顆で求める方法について図１２を用いて説明する。大腿骨-脛骨間の最短距離は外側顆、内側顆それぞれにおいて、骨として抽出された全画素について大腿骨-脛骨間の距離が最小になる点の組合せを全数探索により求める（図１２（ａ）参照）。ただし、このようにして求めた最短距離Ｄｍｉｎの点は量子化による変動を強く受けるため、得られた最短距離から多少の余裕δをみた距離
の範囲に入る組み合わせをすべて用い、それらに距離に応じた重み付けをした上で重心を算出し、最近接点とする（図１２（ｂ））。

以上の処理を、各時系列的に並べられた３次元ＣＴ画像データそれぞれに行なうことにより、最接近点の時間推移を追うことが可能になる（図１２（ｃ））。さらに、上で述べた屈曲角を横軸にとり、最接近点間の距離を縦軸にとったグラフを描くことできる。なお、この場合、最接近点の座標も同時に取得できているため、この点の推移を時系列的に表示して動画像として表示するといったことも可能である。

以上、４次元ＣＴ装置等の４次元撮像装置で得られた骨や軟骨と神経、血管などの軟部組織の３次元連続画像（動態画像）を各要素画像に細分化し、それらの各画像間で同一部分である構成要素をマッチングさせることで構成要素の変位を解析できる。骨などの硬組織では、位置情報の変化（動きの変化）を解析でき、軟部組織ではどのように形状が変化していくかを解析できることで、病気の診断、治療方法の決定、治療効果の判定に使用することができる。

なお、本実施形態では４次元ＣＴ装置を用いた例を説明したが、上述のとおり、同様の手法により画像を解析することができる限りにおいて４次元ＭＲＩ装置など他の撮像装置により得られた４次元画像データに対しても同様に適用が可能である。

４次元ＣＴ装置の原理を示すための図。 人間の膝関節における３次元ＣＴ画像データの模式図。 プログラムの処理を示すフローチャート図 ３次元ＣＴ画像データ構成について説明する図。 開始点探索のアルゴリズムを説明するための模式図。 膝関節における抽出対象領域を分割した結果について示す図。 領域分割の結果から得られる表面画を示す図。 変位量の算出の原理を示す概念図。 軟体に対し非線形変換のパラメータを更に加えることで変位量を得る場合の概念図。 局所的な変形の情報と荷重との関係についての概念図。 膝の屈曲角に対する大腿骨−脛骨間の最近接点間距離および最近接点位置の推移をグラフ化する方法について説明する図。 大腿骨-脛骨間の最短距離を外側顆、内側顆で求める方法について説明する図。

符号の説明

１…Ｘ線源、２…２次元検出器、３…被写体、４…リング、５…開始点

Claims (7)

1. コンピュータに、時系列的に得られた複数の３次元画像データのそれぞれに対して抽出対象領域を分割するステップ、前記分割された複数の前記抽出対象領域における変位量を求めるステップと、を実行させるための画像処理プログラムであって、
   前記時系列的に得られた複数の３次元画像データのそれぞれにおいて抽出対象領域を分割するステップは、前記３次元画像データにおける最大の画素値を求めるステップ、前記求めた最大の画素値よりも小さい所定の画素値を設定するステップ、前記３次元画像データを複数の２次元画像データに分割してそのいずれかを選択し、該選択された２次元画像データのいずれか一側辺から順次走査して前記所定の画素値以上となった画素を開始点と定め、当該開始点となる画素と隣接する画素における画素値が前記所定の画素値以上であると判断した場合、更に当該画素に隣接する画素の画素値に対し同様に判断を繰返すステップ、を有し、
   前記分割された前記抽出対象領域の変位量を求めるステップは、前記時系列的に得られた複数の３次元画像データそれぞれに対し、分割された前記抽出対象領域をスライスする２次元画像データを複数選択し、そのそれぞれの２次元画像データにおける前記抽出対象領域の重心点を求め、求めた前記複数の重心点を用いて主成分分析を行い前記抽出対象領域の軸を定め、当該軸の変位量を求めるステップ、を有することを特徴とする画像処理プログラム。
2. 前記３次元画像データは、ＣＴ画像データである請求項１記載の画像処理プログラム。
3. 前記変位量は、回転及び平行移動の量に基づいて求められることを特徴とする請求項１記載の画像処理プログラム。
4. 前記変位量は、非線形変形のパラメータにも基づいて求められることを特徴とする請求項３記載の画像処理プログラム。
5. 前記抽出対象領域は３次元画像データのそれぞれにおいて複数抽出されることを特徴とする請求項１記載の画像処理プログラム。
6. 前記分割された複数の前記抽出対象領域の変位量を求めるステップは、前記時系列的に得られた複数の３次元画像データそれぞれに対し、分割された前記抽出対象領域をスライスする２次元画像データを複数選択し、そのそれぞれの２次元画像データにおける前記抽出対象領域の重心点を求め、求めた前記複数の重心点から前記抽出対象領域の軸を定め、当該軸の変位量を求めるステップ、を有することを特徴とする請求項５記載の画像処理プログラム。
7. 前記分割された複数の前記抽出対象領域の時系列的な変位量を求めるステップは、更に、前記時系列的に得られた複数の３次元画像データそれぞれにおいて、分割された前記抽出対象領域同士の最短距離を求めるステップ、を有することを特徴とする請求項５記載の画像処理プログラム。