まえがき
放射線医学の一部門である核医学における動態機能解析の方法として,連続する画像*)の各画素の動的変化を多変量解析することにより,いくつかの基本的バターンを自動的に抽出する,画像工学におけるKarhunen-loeve変換と類似の手法があり,ファクター・アナリシス(Factor Analysis)と呼ばれている**)。
連続する部分画像に対してこの解析法を行うという新しい処理法を開発した。これによって因子画像の連続的変化から,呼吸機能としての横隔膜,胸壁などの動きが求められ,それらの画像の相関の時間的変化を求めると,正常例と病的例とでは異なったパターンを示し,臨床診断上の有意な情報が得られた。
内容
Ⅰ.原理とアルゴリズム
従来の動態機能解析の方法は図1に示すように各画素(領域)値の時間的変化(タイム・アクティビティ・カーブ,TAC)を,例えば,コンパートメント・モデルに従うとしてパラメーターを求め,その値をファンクショナル・イメージとして表示するのが一般的である。しかしながら,隣合う画素(領域)に対する考慮がされていない。
これに対し,ファクター・アナリシスは全領域の相関をもとに分析している点で画期的である。この方法の原理は連続する画像にはいくつかの生理学的な時間的動きの因子(qを因子の種類として,φq(t)とする)があり,それがどこにどの程度多く含まれているかを因子画像Fq(x,y)で表すことである。
即ち,x,yを位置,tを時間の指標とする連続画像X(x,y;t)を,Qをファクターの数(解析者が設定する)とし,
と近似して,この両辺間の残差が最小になるような,Fq(x,y)とφq(t)を求めるものである。概念図を図2に示す。
これは各ピクセル問の相関係数を要素とした相関行列の固有値を求め,固有値の大きい順に対応する固有ベクトルを採るという,Karhunen-Loeve(KL)変換と基本的に同じ手法である。基本的な手法はD.C.Barbar1)の方法と同じであるが,Qをファクターの数としたとき,(Q-1)個の固有ベクトルを採り(Q-1)次元空間でそれらを総て取り囲むようなQ点を求めると,これらのQ個の点からの各ベクトルの結合定数が非負であること,及びそれがちょうど分配係数になっていることを利用して因一子画像Fq(x,y)が求まり,次いで最小二乗法によりφq(t)が求められる。
問題は求められた因子画像が常に同じ順序で同じ生理学的因子を示すとは限らず,その解釈が解析者の主観によらなければならず,ファクターの数が多くなる程,また複雑な動きをするような疾患例になる程,困難になることである。そのため自動識別させるのには多くの解析例をデータベース化して人工知能的技法(エキスパート・システム)による方法が必要になるだろうとされていた。
この問題を解くのがこの課題の最大の要点であるが,解析の過程を詳しく検討すると,この解決法の一手法が見通せたので以下に,概略を述べる。画素のタイム・アクティビティ・カーブ(TAC)を求め,全領域にわたるTACの相関係数を要素とする行列の固有値を求め,最大固有値から順に対応する固有ベクトルをとる。ファクターの数をQとしたとき,KL変換ではQ個までとるが,負の成分を回避するため(Q-1)個までとする。
次の行列で表される方程式を解く。
例えば,3因子の場合,図4に示すように2次元平面に(Yl,γ2i)がプロットされ,Xiとする。
において,X、が△STUの内部にあれば
si,ti,uiの全てが正または0
そうでなければsi,ti,ui,のいずれかが負
いずれにしてもsi+ti+ui=1がなりたつ。
そこで,全てのXiを取り囲む△STUを求める。この方法としては,図5に示すように各点問の距離の合計が大きい順に3点を選んで三角形を決め,各辺から外側に最も離れた点まで平行移動して△STUを求めている***)。
この方法が最善かどうか,検討の余地があるが,計算時間の短縮のため採用している。このようにして(β11,β21),……,(β13,β23)を求めると上の方程式が解けて,
なるαが得られる。
つまり,al,α2,,asiは各画素がどの因子に近いかを0から1までで表す,全部足すと1になる規格化された分配係数になっていて,しかも画像として表せる非負の値になっている。
ここで自動識別する処理が働く。即ち,各頂点から最も距離的に近いX、を求め,それがどの画素かをみる。そうするとその頂点からのαが最も大きくなっているから,その頂点を示すインデックスをその画素(解剖,生理学的な因子に対応する)のものに置き換える。図5に示すように,例えば,横隔膜,胸壁,死腔の因子の順に求めたいのであれば,対応する画素から最短距離の頂点を示すβのインデックスを置き換えれば,よいことになる。
4因子以上の場合も同様にして,(因子の数一1)元空間のプロットを取り囲む(因子の数)個の頂点をもつ超立体を求め,頂点のインデックスを置き換えることによって得られる。
元画像列の平均値に,こうして得られたαの値(の分布)を掛けると因子画像Fq(x,y)が求められる。次いで最小二乗法によりφq(t)が求められる。
II.データ収集と解析法
81mkr一ガス肺換気シンチグラムを呼吸同調装置3)を用いてリスト・モードで収集し,心電同期の場合と同じ乎法で,呼期開始から始まり吸期末端で終わる16枚の周期的連続画像に変換した。図3に示すように,連続画像全体ではなく,1枚目から3枚目,次に2枚目から4枚目,……,の様に部分的に連続する3枚の画像に対して,ファクター・アナリシスを行い,因子画像を自動識別して繋ぐ方法(仮称Partial Sequence Factor Analysis PSFA)を試みた。
周期的連続画像であることから,16枚目と1枚目,2枚目の画像に対する因子画像を1番目とし,1～3枚目に対するものを2番目,以下同様にして15枚目と16枚目,1枚目に対するものを16番目の因子画像として16枚の周期的連続因子画像を求めた。
臨床画像処理装置で解析法を開発,実行するには,時間と場所の問題から日常の検査業務に差し障りがあるので,磁気テープを介して,パソコンにデータを転送し,解析法の開発,実行を続けた。開発言語はCを用いた。
画素数は64×64とした。ただし8×8画素を1つにまとめ相関行列は64×64行列とした。
成果
従来のファクター・アナリシスの評価法は各因子画像のカウント数が全因子画像のカウント総数に対してどの位の割合で寄与しているかであった。しかしながら,図6に示すように,同じ寄与率を示しても因子画像の分れ方が異なる場合が有りうる。したがって,この評価法よりも,呼吸機能としての横隔膜と胸壁との間でいかにアイソトープがトレーサーとして流れているかを二つの因子画像がいかに分離しているか,又は重なっているかという点から両者の画像間の相関係数で評価した。
このようなトレーサーとしてのアイソトープの流れが一呼吸の間でどのようになっているのかを見いだすためPSFAによる解析を行った。一呼吸の間の横隔膜と胸壁の動きの正常者の例を図8に,閉塞性肺疾患例の例を図9に示す。いずれも後面座位像である。
正常例では,呼期の途中の4番目から6番目にかけてと,吸期の途中の12番目から14番目にかけて横隔膜像と胸壁像とがよく分離されているが,この時期を過ぎると2つの像の分離は悪くなる。
各時点での因子画像間の相関係数の変化を求め,図7に示す。正常例では吸期から吸期に,逆に呼期から呼期に転ずる時点で相関が高くなっており,ここで横隔膜と胸壁との間のガス交換が行われていると考えられる。
一方,閉塞性肺疾患例の例では一呼吸期間全体にわたり相関が高くなっていて,正常例のようなガス交換ができていないことが示唆される。
因子画像は同一被検者の,同一部位の同一時点のものであるから,解剖学的位置関係がよくわかるように,2ファクターの場合の二つの因子画像の重ね合わせ表示(通常は行われないシャウカステン上で,二枚のフィルムを重ね合わせて見ることを,CRTディスプレイ.上で実現すること)を試みた。図10に,図8,9の因一子画像の動きを重ね合わせた例を示す。因一子画像は元の画像を規格化された分配係数の分布によって分配したものであるから単純に加算すると元画像になる。そこで,1画素毎に交互に二つの因子画像を表示した。見やすいようにカラースケールは変えている。
正常例では,息を吐き始める時点では肺の外側がよく動き,逆に息を吸い始める時点では内側(気管支付近)がよく動いていることが得られたJ一方,閉塞性肺疾患例の例では横隔膜と胸壁の区別をつけ難い,複雑な動き
を示す。更に解析を細分化していくと,呼吸機能の生理学・病理学な新たな情報,例えば「閉塞がどの部位,時点で生じているのか」などが得られると思われる。
まとめ
この手法の考え方を簡潔に表現すると,「画像圧縮に用いられる,統計的直交変換の変換核となる基底を画像として求め,その動きをとらえる」ということになる。
元画像列を,単に前後に差分を取るなどの手法では,このような因子画像の動きは得られない。局所的な因子画像列を繋ぐことは,臨床研究者にとって興味深い「広範囲に分布するRIの中で,ある時間のRIが次の瞬間にどこへ移動したか」を求める一つのアプローチになると思われる。
重ね合わせ表示については,さらに見やすい(因子画像の区別がつきやすい)ように改良が望まれる。