逐次近似法の基礎と各種補正方法

逐次近似法の基礎と各種補正方法 横浜創英大学橋本雄幸 画像再構成における逐次近似法の歴史は長く,X 線 CT においても解析的方法が見つかる前は, 逐次近似法を用いて画像を再構成していた. 解析的方法が見つかってからは, 計算時間の長さから逐次近似法はあまり使われなくなった. しかし, コンピュータの発展に伴い, 繰り返しても計算時間がそれほどかからなくなったこともあり, 解析的方法が確立できない SPECT の画像再構成や線量を抑えた X 線 CT においても逐次近似法が再び用いられるようになった. ここでは, その逐次近似法の基本的な考え方から出発し, 画像再構成に使われている逐次近似法の意味まで紐解いていく. 1. 逐次近似法の考え方問題を解く一般的な手順は, 計算の道具である数式を実際の計測に合わせて定式化するところから始まる. 実際の計測を定式化することを 順問題 という. 次に, 順問題によって定式化した式を解くことによって答えを導く. 順問題の式から逆に解くことを 逆問題 という. 例えば, 10 個のリンゴを 0 g の皿に載せて重さを測ったら 10 g になった. という計測があったとする. ここからリンゴ 1 個あたりの重さを とし, 測った重さを として定式化するのが順問題である. その式は 10 個のリンゴと 0 g の皿 から 10 0 (1) となる. これを について解析的に解くのが逆問題で,(1) 式を について解くと 0 10 となる. この 順問題 と 逆問題 の関係を図 1 に示す. 式の に 10 を代入すれば, リンゴ 1 個あたりの重さが 100 g であることが直接求められる. 式のように逆問題を解析的に解くことによって,10 個のリンゴの計測結果が分かれば, リンゴ 1 個あたりの重さはすぐに求めることができる. では,(1) 式の順問題だけで解を求めることを考えよう. に適当な数字を当てはめて順問題の式が成り立てば, そのときの が解である. の決め方はいろいろ考えられるが, 最も単純な方法は, 当てはめる の値を小さな値から順に大きくしていき, 順問題の結果が等しくなるかまたは最も近くなるところを解とすることができる. 例えば, 最初に 10 を当てはめると順問題の結果は 1010 0 300 となる. さらに,30, と当てはめていくと

10 0 400 10 30 0 500 10 40 0 600 10 50 0 700 10 60 0 800 10 70 0 900 10 80 0 1000 10 90 0 1100 10 100 0 10 (4) となり,100 を当てはめたところで結果が一致する. このように解を求めることもできる. 次にもう少し効率的な求め方を考える. 例えば 式のように に 10 を当てはめたとき, 順問題の結果は 300 となり, 計測値の 10 とは 900 の差が出る. その差を利用して次の を決めることを考える. これは制御でよく用いられるフィードバックの考え方である. フィードバックは図 2 に示すような形で表すことができる. 最初の値 (0) = 10 とすると計測値 との差は (0) (10 0) 10 (1010 0) 10 300 900 (5) と求められる. 次の を決めるのに差の 900 を利用するが, そのまま最初の値に加えて次の値を決める と徐々に値が発散してしまうので, 戻す値を小さく調整する必要がある. 例えば 1/ にして戻すことを 考えると,900 / = 45 を戻すことになるので, 次の値 (1) は (0) (1) (0) (10 0) 900 10 10 45 55 (6) となる. これを繰り返すと, (4) (5) (1) (4) (1) (10 0) 10 (10 55 0) 55 77.5 (10 0) 10 (10 77.5 0) 77.5 88.75 (7) (10 0) 10 (10 88.75 0) 88.75 94.375 (4) (10 0) 10 (10 94.375 0) 94.375 97.1875 となり, 徐々に解の 100 に近づくことが分かる. ここで,k 番目の の値 (k) から k+1 番目の (k+1) を求め るときのフィードバックの値は, 図 2 で示しているように

( k (10 ) 0) (8) となるので, 一般式は ( k1) (10 0) (9) と表現できる. このように順問題だけで繰り返しながら解に近づける方法が逐次近似法である. ここで示したのは加減型の逐次近似法で, 計測値と順問題の結果との差分を出して, 仮定した値に加えて次の値を出している. この模式図を図 3 に示す. 加減法では, 値を戻すときにその値を調整する必要があり, 差分に調整値を乗算する.(9) 式では,1/ を調整値として乗算していることになる. 調整値を 1/ にしたときの の値の収束の様子を表したグラフを図 4 に示す. 調整値をαとしたときの一般式は ( k1) ( k) ( k) [ (10 0)] (10) となる. この調整値を 1/6,1/5,1/4 にしたときの の値の変化を表したグラフをそれぞれ図 5~7 に示す. 図 5 の調整値が 1/6 の場合は値が上下しながらも収束している. しかし, 図 6 の調整値が 1/5 の場合は値が上下に振動してしまい, 図 7 の調整値が 1/4 の場合は値が上下しながら発散してしまう. 調整値 αは小さい程とゆっくり収束するが, 大きいと変化が激しくなり, 値によっては収束せずに発散してしまう. このように加減型の場合は, 調整値によって収束の様子が変化する. 逐次近似法には, 図 8 に示したような乗除型のフィードバックを利用したものもある. 乗除型の場合, フィードバックの値は計測値と順問題の値の比 ( 除算 ) となるので, 式で表すと 10 ( k ) 0 (11) となり, 一般式は ( k1) (12) 10 0 と表現できる. 乗除型の模式図を図 9 に示す. 乗除型では調整値の必要がなく, 比と乗算によって次の 値が求まる. 初期値 (0) を 10 とすると,

(1) (4) (0) (1) 10 10 10 10 10 10 40 0 10 10 0 10 40 80 0 10 40 0 10 80 96 0 10 80 0 10 96 99.31 0 10 96 0 (0) (1) (13) となり, この場合も解の 100 に近づくことが分かる. (12) 式をもとにした の値の収束の様子を図 10 に示す. 乗除型の場合は調整値がなくても安定して収束する傾向がある. 2.Ecel を利用した逐次近似法のシミュレーション Ecel を利用して逐次近似法の収束の様子をシミュレーションする. 順問題には前節の 10 個のリンゴを 0 g の皿に載せて重さを測ったら 10 g になった. を利用する. この節では, 前節で示した収束のグラフを作成することを目標とする. まずは加減型の逐次近似法についてシミュレーションを行なう. 図 11 に示すようにあらかじめ Ecel に文字列と値を入力する. 調整値のセル B1 には, 調整値の逆数である分母の値 を入力する. 繰り返し回数の 0 回目のセル B3 は初期値になるので, 前節で設定した初期値の 10 を入力しておく. セル B4 には(9) 式に対応する式を入力する.(9) 式において (k) はセル B3 となり, は計測値の 10, 分母の はセル B1 となる. セル B1 は複写の際に参照セルをずらさないようにするので, 絶対参照の $B$1 とする. よって, セル B4 には図 12 に示すように =B3+(10-(B3*10+0))/$B$1 (14) と入力する. ここで,(9) 式と (14) 式が対応していることを確認していただきたい. Enter キーを押して決定すると, 55 という計算結果が表示される. このセル B4 を下方向に複写すと, 図 13 に示すような数値が計算される. この結果をもとに散布図を作成する手順を以下に示す. この手順は Ecel10 をもとにしている. 1 繰り返し回数と の値のセル範囲 A3:B18 を選択する. 2 挿入 タブの 散布図 散布図( 直線とマーカー ) をクリックする. 図 4 のような散布図のグラフが表示される. セル B1 の調整値を変えることによって, 図 5~7 のような収束の変化を観察できる. また, 調整値を様々に変えて, その変化を観察して欲しい. 次に乗除型の逐次近似法についてシミュレーションを行なう. 図 11 と同様にあらかじめ Ecel に文字列と値を入力する. ただし, 調整値は必要ないので省いておく. 繰り返し回数の 0 回目のセル B3 には先ほどと同様に初期値 10 を入力する. セル B4 には(12) 式に対応する式を入力する.(12) 式において (k) はセル B3 となり, は計測値の 10 となる. よって, セル B4 には図 14 に示すように

=B3*10/(B3*10+0) (15) と入力する. ここでも,(12) 式と (15) 式が対応していることを確認していただきたい. Enter キーを押 して決定すると, 40 という計算結果が表示される. このセル B4 を下方向に複写すと, 図 15 に示 すような数値が計算される. 加減型の場合と同様に, この結果をもとに散布図を作成することができる. 3. 画像再構成における逐次近似法画像再構成における逐次近似法は ART 法や SIRT 法と呼ばれる加減型が最初に提案された.X 線 CT の画像再構成における順問題と逆問題は図 16 に示すように, 順問題の一般式はラドン変換となり, 逆問題を解析的に解いたものがフーリエ変換法や FBP 法となる. 逐次近似法は順問題の式のみで解いていく方法なので, 使う式は投影に当たるラドン変換となる. 被写体の画像を f(,, ラドン変換によって投影したデータを g(x,θ) とすると, ラドン変換の式は g ( X, ) dy (16) となる. ラドン変換は Y 軸の方向に積分する形になるが, これは被写体の画像をその方向に足していくことに相当する. ラドン変換によって投影することで, 画像空間 (, から投影空間 (X,θ) へと空間が変換される. フィードバックを利用するときに, 投影空間から画像空間へ空間を戻す必要があり, 逆投影という手法を利用する. 逆投影の式は, 1 2 g( X, ) d (17) 2 0 となる. 逆投影は角度方向に積分する形になるが, 投影データをあらゆる角度方向から足し合わせるこ とに相当する. 実際には 1 周の 2π で割っているので, 平均を出していることになる. 投影と逆投影の 関係を図 17 に示す. 逆投影によって画像空間には戻せるが, 逆投影は全角度方向への平均値を求める 計算なので, 逆投影で求まった画像は非常にぼけたものになる. 加算型の逐次近似法をこの画像再構成に当てはめると次のようになる. 初期値 : (0) 初期画像 f(, (0) ( すべて 1 の画像など ) 順問題:10 (k) +0 計測値 : g(x,θ) (k) dy ( k) フィードバック: [ (10 0)] ( k) [ g( X, ) dy] の逆投影 ここで α はどちらも調整値に相当する. 画像再構成のフィードバックには空間を戻すための逆投影が 含まれるが, 基本の形は変わらない. フィードバックを図に表したものを図 18 に示す. これを k 番目

の繰り返しに当てはめて再構築すると f 1 2 2 ( k1) (, [ g( X, ) dy] d 0 (18) となる. この式は SIRT 法と呼ばれる逐次近似法に相当する. 実際に数値ファントムの画像を利用して 計算機シミュレーションを行なった結果を図 19 に示す. ここで, 調整値 α は 1/1000 を用いている. 加 算型の逐次近似法では, 画像再構成においても, 調整値 α の値によって収束の様子は変化する. 調整値 α が 1/900 の場合を図 に示す. この場合は値が発散して, ある程度繰り返すと画像が乱れてしまう. 次に, 乗除型の逐次近似法を画像再構成に当てはめると次のようになる. 初期値 : (0) 初期画像 f(, (0) ( すべて 1 の画像など ) 順問題:10 (k) +0 計測値 : g(x,θ) (k) dy フィードバック: 10 ( k ) 0 g( X, ) (k ) dy の逆投影 ここでも画像再構成のフィードバックには空間を戻すための逆投影が含まれるが, 基本の形は変わら ない. フィードバックを図に表したものを図 21 に示す. これを k 番目の繰り返しに当てはめて再構築 すると f g( X, ) 2 ( k1) (, 2 d 0 dy 1 (19) となる. 実際に数値ファントムの画像を利用して計算機シミュレーションを行なった結果を図 22 に示す. また, 加減型のα=00,1000,900 の場合と乗除型の逐次近似法で再構成した画像の評価値を, 繰り返し回数でグラフにした結果を図 23 に示す. 加減型は調整値 αの値によって収束の様子が異なるのに対し, 乗除型は調整値がなくても安定していることが分かる. 4.ML-EM 法による画像再構成 SPECT の画像再構成などに使われている逐次近似法の 1 つである ML-EM 法では, 検出確率 C i という概念を導入する. まずは,SPECT 画像再構成における減弱, 散乱, コリメータ特性などの問題がない場合を仮定する. この場合,X 線 CT の画像再構成法と等価になる. 検出確率 C i は, 番目の画素から放出された放射線が i 番目の検出器で検出される確率を意味する. 検出器の前にはコリメータがあり, 理想的な場合は検出器に垂直に入った放射線のみがカウントされる. 減弱などがなければ, 単純な投影と同じことになる. ある画素から放出された放射線は特定の検出器にしか入らないので, 検出確率を画像で見ると図 25 に示すようにほとんどが 0 となり, 値が入るところは細かい直線が並んだようになる.

検出確率 C i を用いた順問題 ( 投影 ) の式は M 1 i C i 0 () と表される. ここで λ は画像で i は投影データである. この式は,(16) 式のラドン変換と等価であり, 図 26 に示すような行列と列ベクトルとの掛け算になる. また, 検出確率 C i を用いた逆投影の式は N 1 1 ' 1 ic N i (21) i0 C i0 i と表される.() 式の投影とは C i を逆に掛けた形になっている. この式は,(17) 式の逆投影と等価であり, 図 27 に示すような行ベクトルと行列との掛け算になる. ここで分数は (17) 式の 1/2πと同様, 平均値を出すための演算である. ML-EM 法では k 番目の画像をλ (k),k+1 番目の画像をλ (k+1) とするとき, 以下の式が用いられる. N 1 ( k1) 1 i Ci (22) N 1 C C i0 i M 1 i0 ' 0 i' ' この式の導出は, 他の書籍を見ていただきたい. ここでは, この式の意味するところを考える. フィー ドバックの形で表すと図 28 に示すようになる. この形は図 21 の乗除型と等価であることが分かる. よ って,ML-EM 法は乗除型の逐次近似法である. その対応を細かく見ると以下のようになる. 初期画像 ( すべて 1 の画像など ):f(, (0) λ (0) 順問題( 投影 ): 計測値 :g(x,θ) i フィードバック : M 1 (k) dy C g( X, ) (k ) 0 の逆投影 dy i M 1 ' 0 C i i' ' の逆投影 実際に減弱などのない ML-EM 法で行った結果は,X 線 CT の乗除型の逐次近似法と同様の結果となる. ここで重要になるのが検出確率 C i の存在である. 実際の SPECT の順問題は, 図 29 に示すように, 投影の他に減弱や散乱, コリメータ特性などの問題が入ってくる. その逆問題は難解となるので, でき れば順問題を確実に定式化し, 順問題のみで解ける逐次近似法を用いることが望まれる.SPECT の計測 過程は図 30 に示すように, 検出確率 C i を用いることでうまく定式化することが可能である.C i を求め れば (22) 式の繰り返しの式によって解である画像を求めることができる. 実際には C i にどこまで ( 減弱,

散乱, コリメータ特性 ) の問題を含めるかや C i の求め方などは, メーカーや装置によって異なるが, 逐次近似法の考え方は変わらない. 最近,X 線 CT にも線量を少なくするなどの目的で,ML-EM 法を画像再構成に用いる例が見られるようになった.X 線 CT などの透過型 CT の ML-EM 法は N 1 M 1 li ( di ep li i ) ( k 1) i0 0 (23) N 1 M 1 M 1 li li di ep li i0 0 0 と表され, 図 31 にも示すように, 加減型の逐次近似法となる. 加減型においてはフィードバックの値を調整しなければならなく,(23) 式においてもωなどの調整値が存在する. よって, 画像を収束させるためには,ωなどの調整値をうまく設定する必要がある. 以上のように, 繰り返しの方法がどのようなことを行っているのかを, その意味合いから知ることは重要である. 実際の装置では, 逐次近似法の基本式に加え, 様々な工夫が盛り込まれている. その工夫の部分を紐解くことはできないが, 逐次近似法の基本をつかんでいただければ, 少しでもブラックボックスが解消されることでしょう. 特に加減型と乗除型の性質の違いだけでも感じてもらえれば幸いです. 図の説明 図 1 順問題と逆問題の関係 図 2 加減型のフィードバックの考え方

図 3 加減型の逐次近似法の模式図 図 4 差分に 1/ を掛けてから加える という フィードバックを行って繰り返した結果 図 5 差分に 1/6 を掛けてから加える という フィードバックを行って繰り返した結果 図 6 差分に 1/5 を掛けてから加える という フィードバックを行って繰り返した結果 図 7 差分に 1/4 を掛けてから加える という フィードバックを行って繰り返した結果 図 8 乗除型のフィードバックの考え方

図 9 乗除型の逐次近似法の模式図 図 10 比をとってから掛ける という フィードバックを行って繰り返した結果 図 11 事前に入力しておく文字列と値 図 12 加減型の繰り返しの式の入力 図 13 加減型の繰り返しの計算結果 図 14 乗除型の繰り返しの式の入力

図 15 乗除型の繰り返しの計算結果 図 16 X 線 CT の順問題と逆問題 図 17 投影と逆投影の関係 図 18 X 線 CT の加算型の逐次近似法を フィードバックで示した図 図 19 加算型の逐次近似法で数値シミュレ ーションを行なった結果 (α=1/1000) 図 加算型の逐次近似法で α=1/900 にして数値 シミュレーションを行なった結果

図 21 X 線 CT の乗除型の逐次近似法を フィードバックで示した図 図 22 乗除型の逐次近似法で数値シミュ レーションを行なった結果 図 23 加減型と乗除型の繰り返し回数と 再構成画像の評価値 図 24 検出確率 C i の考え方 図 25 16 16 画素の画像と 16 16 投影の 図 26 検出確率 C i と投影の計算 投影データから作成した検出確率 C i の画像 (256 256 画素 )

図 27 検出確率 C i と逆投影の計算 図 28 ML-EM 法をフィードバックで表した図 図 29 SPECT の順問題と逆問題 図 30 検出確率 C i による定式化 図 31 透過型 CT の ML-EM 法をフィードバック で表した図 ( 加減型となる )