逐次近似法の基礎と各種補正方法
横浜創英大学 橋本雄幸
画像再構成における逐次近似法の歴史は長く，X線CTにおいても解析的方法が見つかる前は，逐次 近似法を用いて画像を再構成していた．解析的方法が見つかってからは，計算時間の長さから逐次近似 法はあまり使われなくなった．しかし，コンピュータの発展に伴い，繰り返しても計算時間がそれほど かからなくなったこともあり，解析的方法が確立できないSPECTの画像再構成や線量を抑えたX線CT においても逐次近似法が再び用いられるようになった．ここでは，その逐次近似法の基本的な考え方か ら出発し，画像再構成に使われている逐次近似法の意味まで紐解いていく．
1．逐次近似法の考え方
問題を解く一般的な手順は，計算の道具である数式を実際の計測に合わせて定式化するところから始 まる．実際の計測を定式化することを「順問題」という．次に，順問題によって定式化した式を解くこ とによって答えを導く．順問題の式から逆に解くことを「逆問題」という．
例えば，「10個のリンゴを200 gの皿に載せて重さを測ったら1200 gになった．」という計測があった とする．ここからリンゴ1個あたりの重さをxとし，測った重さをyとして定式化するのが順問題であ る．その式は「10個のリンゴと200 gの皿」から
200 10
x
y
(1)となる．これをxについて解析的に解くのが逆問題で，(1)式をxについて解くと
10
200
y
x
(2)となる．この「順問題」と「逆問題」の関係を図1に示す．(2)式のyに1200を代入すれば，リンゴ1 個あたりの重さが100 gであることが直接求められる．(2)式のように逆問題を解析的に解くことによっ て，10個のリンゴの計測結果が分かれば，リンゴ1個あたりの重さはすぐに求めることができる．
では，(1)式の順問題だけで解を求めることを考えよう．xに適当な数字を当てはめて順問題の式が成 り立てば，そのときのxが解である．xの決め方はいろいろ考えられるが，最も単純な方法は，当ては めるxの値を小さな値から順に大きくしていき，順問題の結果が等しくなるかまたは最も近くなるとこ ろを解とすることができる．例えば，最初に10を当てはめると順問題の結果は
300 200 10
10
y
(3)となる．さらに20，30，・・・と当てはめていくと
1200 200
100 10
1100 200
90 10
1000 200
80 10
900 200 70 10
800 200 60 10
700 200 50 10
600 200 40 10
500 200 30 10
400 200 20 10
y y y y y y y y y
(4)
となり，100を当てはめたところで結果が一致する．このように解を求めることもできる．
次にもう少し効率的な求め方を考える．例えば(3)式のようにxに10を当てはめたとき，順問題の結 果は300となり，計測値の1200とは900の差が出る．その差を利用して次のxを決めることを考える．
これは制御でよく用いられるフィードバックの考え方である．フィードバックは図2に示すような形で 表すことができる．最初の値x(0) = 10とすると計測値yとの差は
900 300 1200 )
200 10 10 ( 1200 )
200 10
(
(0)
x
y
(5)と求められる．次のxを決めるのに差の900を利用するが，そのまま最初の値に加えて次の値を決める と徐々に値が発散してしまうので，戻す値を小さく調整する必要がある．例えば1/20にして戻すことを 考えると，900 / 20 = 45を戻すことになるので，次の値x(1)は
55 45 20 10
10 900 20
) 200 10
(
(0)) 0 ( ) 1
(
y x
x
x
(6)となる．これを繰り返すと，
1875 . 20 97
) 200 375 . 94 10 ( 375 1200
. 20 94
) 200 10
(
375 . 20 94
) 200 75 . 88 10 ( 75 1200
. 20 88
) 200 10
(
75 . 20 88
) 200 5 . 77 10 ( 5 1200 . 20 77
) 200 10
(
5 . 20 77
) 200 55 10 ( 55 1200
20
) 200 10
(
) 4 ( )
4 ( ) 5 (
) 3 ( )
3 ( ) 4 (
) 2 ( )
2 ( ) 3 (
) 1 ( )
1 ( ) 2 (
x x y
x
x x y
x
x x y
x
x x y
x
(7)
となり，徐々に解の100に近づくことが分かる．ここで，k番目のxの値x(k)からk+1番目のx(k+1)を求め るときのフィードバックの値は，図2で示しているように
20
) 200 10
(
( )
x
ky
(8)となるので，一般式は
20
) 200 10
(
( )) ( ) 1
(k
k y x
k
x
x
(9)と表現できる．このように順問題だけで繰り返しながら解に近づける方法が逐次近似法である．ここで 示したのは加減型の逐次近似法で，計測値と順問題の結果との差分を出して，仮定した値に加えて次の 値を出している．この模式図を図3に示す．加減法では，値を戻すときにその値を調整する必要があり，
差分に調整値を乗算する．(9)式では，1/20 を調整値として乗算していることになる．調整値を1/20 に したときのxの値の収束の様子を表したグラフを図4に示す．
調整値をαとしたときの一般式は
)]
200 10
(
[
( )) ( ) 1
(k
x
k y x
k
x
(10)となる．この調整値を1/6，1/5，1/4にしたときのxの値の変化を表したグラフをそれぞれ図5～7に示 す．図5の調整値が1/6の場合は値が上下しながらも収束している．しかし，図6の調整値が1/5の場 合は値が上下に振動してしまい，図7の調整値が1/4の場合は値が上下しながら発散してしまう．調整 値αは小さい程とゆっくり収束するが，大きいと変化が激しくなり，値によっては収束せずに発散して しまう．このように加減型の場合は，調整値によって収束の様子が変化する．
逐次近似法には，図8に示したような乗除型のフィードバックを利用したものもある．乗除型の場合，
フィードバックの値は計測値と順問題の値の比（除算）となるので，式で表すと
200 10 x
(k)
y
(11)となり，一般式は
200 10
( )) ( ) 1 (
k k
k
x x y
x
(12)と表現できる．乗除型の模式図を図9に示す．乗除型では調整値の必要がなく，比と乗算によって次の 値が求まる．初期値x(0)を10とすると，
31 . 200 99 96 10 96 1200 200
10
200 96 80 10 80 1200 200
10
200 80 40 10 40 1200 200
10
200 40 10 10 10 1200 200
10
) 3 ( )
3 ( ) 4 (
) 2 ( )
2 ( ) 3 (
) 1 ( )
1 ( ) 2 (
) 0 ( )
0 ( ) 1 (
x x y
x
x x y
x
x x y
x
x x y
x
(13)
となり，この場合も解の100に近づくことが分かる． (12)式をもとにしたx の値の収束の様子を図10 に示す．乗除型の場合は調整値がなくても安定して収束する傾向がある．
2．Excelを利用した逐次近似法のシミュレーション
Excelを利用して逐次近似法の収束の様子をシミュレーションする．順問題には前節の「10個のリン
ゴを200 gの皿に載せて重さを測ったら1200 gになった．」を利用する．この節では，前節で示した収
束のグラフを作成することを目標とする．
まずは加減型の逐次近似法についてシミュレーションを行なう．図11に示すようにあらかじめExcel に文字列と値を入力する．調整値のセル【B1】には，調整値の逆数である分母の値「20」を入力する．
繰り返し回数の0回目のセル【B3】は初期値になるので，前節で設定した初期値の「10」を入力してお く．セル【B4】には(9)式に対応する式を入力する．(9)式においてx(k)はセル【B3】となり，yは計測値
の1200，分母の20はセル【B1】となる．セル【B1】は複写の際に参照セルをずらさないようにするの
で，絶対参照の「$B$1」とする．よって，セル【B4】には図12に示すように
=B3+(1200-(B3*10+200))/$B$1 (14)
と入力する．ここで，(9)式と(14)式が対応していることを確認していただきたい．「Enter」キーを押し て決定すると，「55」という計算結果が表示される．このセル【B4】を下方向に複写すと，図13に示す ような数値が計算される．
この結果をもとに散布図を作成する手順を以下に示す．この手順はExcel2010をもとにしている．
① 繰り返し回数とxの値のセル範囲【A3:B18】を選択する．
② 「挿入」タブの「散布図」→「散布図（直線とマーカー）」をクリックする．
図4のような散布図のグラフが表示される．セル【B1】の調整値を変えることによって，図5～7のよ うな収束の変化を観察できる．また，調整値を様々に変えて，その変化を観察して欲しい．
次に乗除型の逐次近似法についてシミュレーションを行なう．図11と同様にあらかじめExcelに文字 列と値を入力する．ただし，調整値は必要ないので省いておく．繰り返し回数の0回目のセル【B3】に は先ほどと同様に初期値「10」を入力する．セル【B4】には(12)式に対応する式を入力する．(12)式に おいてx(k)はセル【B3】となり，yは計測値の1200となる．よって，セル【B4】には図 14に示すよう に
=B3*1200/(B3*10+200) (15)
と入力する．ここでも，(12)式と(15)式が対応していることを確認していただきたい．「Enter」キーを押 して決定すると，「40」という計算結果が表示される．このセル【B4】を下方向に複写すと，図15に示 すような数値が計算される．加減型の場合と同様に，この結果をもとに散布図を作成することができる．
3．画像再構成における逐次近似法
画像再構成における逐次近似法はART法やSIRT法と呼ばれる加減型が最初に提案された．X線CT の画像再構成における順問題と逆問題は図16に示すように，順問題の一般式はラドン変換となり，逆 問題を解析的に解いたものがフーリエ変換法やFBP法となる．逐次近似法は順問題の式のみで解いてい く方法なので，使う式は投影に当たるラドン変換となる．被写体の画像をf(x,y)，ラドン変換によって投 影したデータをg(X,θ)とすると，ラドン変換の式は
f x y dY X
g ( , ) ( , )
(16)となる．ラドン変換はY軸の方向に積分する形になるが，これは被写体の画像をその方向に足していく ことに相当する．ラドン変換によって投影することで，画像空間(x, y)から投影空間(X,θ)へと空間が変 換される．フィードバックを利用するときに，投影空間から画像空間へ空間を戻す必要があり，逆投影 という手法を利用する．逆投影の式は，
2
0
( , )
2 ) 1 ,
( x y g X d
f
(17)となる．逆投影は角度方向に積分する形になるが，投影データをあらゆる角度方向から足し合わせるこ とに相当する．実際には1周の2πで割っているので，平均を出していることになる．投影と逆投影の 関係を図 17 に示す．逆投影によって画像空間には戻せるが，逆投影は全角度方向への平均値を求める 計算なので，逆投影で求まった画像は非常にぼけたものになる．
加算型の逐次近似法をこの画像再構成に当てはめると次のようになる．
・初期値：x(0) → 初期画像f(x,y)(0)（すべて1の画像など）
・順問題：10x(k)+200 →
f ( x , y )
(k)dY
・計測値：y → g(X,θ)・フィードバック：
[ y ( 10 x
(k) 200 )]
→ [ g ( X , ) f ( x , y )
(k)dY ]
の逆投影
ここでαはどちらも調整値に相当する．画像再構成のフィードバックには空間を戻すための逆投影が 含まれるが，基本の形は変わらない．フィードバックを図に表したものを図 18に示す．これを k 番目
の繰り返しに当てはめて再構築すると
2 0
) ( )
( )
1
(
[ ( , ) ( , ) ]
2 ) 1 , ( )
,
( x y f x y g X f x y dY d
f
k k k (18)となる．この式は SIRT 法と呼ばれる逐次近似法に相当する．実際に数値ファントムの画像を利用して 計算機シミュレーションを行なった結果を図19に示す．ここで，調整値αは1/1000を用いている．加 算型の逐次近似法では，画像再構成においても，調整値αの値によって収束の様子は変化する．調整値
αが1/900の場合を図20に示す．この場合は値が発散して，ある程度繰り返すと画像が乱れてしまう．
次に，乗除型の逐次近似法を画像再構成に当てはめると次のようになる．
・初期値：x(0) → 初期画像f(x,y)(0)（すべて1の画像など）
・順問題：10x(k)+200 →
f ( x , y )
(k)dY
・計測値：y → g(X,θ)・フィードバック：
200 10 x
(k)
y
→
f x y dY X
g
k)
)
(, (
) , (
の逆投影
ここでも画像再構成のフィードバックには空間を戻すための逆投影が含まれるが，基本の形は変わら ない．フィードバックを図に表したものを図 21に示す．これを k番目の繰り返しに当てはめて再構築 すると
2
0 ( )
) ( )
1 (
) , (
) , ( 2
) 1 , ( )
,
( d
dY y x f
X y g
x f y
x f
k k
k (19)
となる．実際に数値ファントムの画像を利用して計算機シミュレーションを行なった結果を図 22 に示 す．また，加減型のα=2000，1000，900の場合と乗除型の逐次近似法で再構成した画像の評価値を，繰 り返し回数でグラフにした結果を図 23 に示す．加減型は調整値αの値によって収束の様子が異なるの に対し，乗除型は調整値がなくても安定していることが分かる．
4．ML-EM法による画像再構成
SPECTの画像再構成などに使われている逐次近似法の1つであるML-EM法では，検出確率Cijとい
う概念を導入する．まずは，SPECT画像再構成における減弱，散乱，コリメータ特性などの問題がない 場合を仮定する．この場合，X線CTの画像再構成法と等価になる．検出確率Cijは，j番目の画素から 放出された放射線がi番目の検出器で検出される確率を意味する．検出器の前にはコリメータがあり，
理想的な場合は検出器に垂直に入った放射線のみがカウントされる．減弱などがなければ，単純な投影 と同じことになる．ある画素から放出された放射線は特定の検出器にしか入らないので，検出確率を画 像で見ると図25に示すようにほとんどが0となり，値が入るところは細かい直線が並んだようになる．
検出確率Cijを用いた順問題（投影）の式は
10 M
j
j ij
i
C
y
(20)と表される．ここでλjは画像で yiは投影データである．この式は，(16)式のラドン変換と等価であり，
図26に示すような行列と列ベクトルとの掛け算になる．
また，検出確率Cijを用いた逆投影の式は
10 1
0
' 1
N
i ij N i
i ij
j
y C
C
(21)と表される．(20)式の投影とはCijを逆に掛けた形になっている．この式は，(17)式の逆投影と等価であ り，図27に示すような行ベクトルと行列との掛け算になる．ここで分数は(17)式の1/2πと同様，平均 値を出すための演算である．
ML-EM法ではk番目の画像をλj(k)，k+1番目の画像をλj(k+1)とするとき，以下の式が用いられる．
10 1
0 '
) ( ' ' 1
0 ) ( ) 1
(
1
Ni M
j
k j ij
ij i N
i ij k
j k
j
C C y
C
(22)この式の導出は，他の書籍を見ていただきたい．ここでは，この式の意味するところを考える．フィー ドバックの形で表すと図28に示すようになる．この形は図21の乗除型と等価であることが分かる．よ って，ML-EM法は乗除型の逐次近似法である．その対応を細かく見ると以下のようになる．
・初期画像（すべて1の画像など）：f(x,y)(0) → λj(0) ・順問題（投影）：
f ( x , y )
(k)dY
→
1
0 M
j j
C
ij
・計測値：g(X,θ) → yi
・フィードバック：
f x y dY X
g
k)
)
(, (
) , (
の逆投影 →
1
0 '
) ( ' ' M
j
k j ij
i
C y
の逆投影
実際に減弱などのないML-EM法で行った結果は，X線CTの乗除型の逐次近似法と同様の結果となる．
ここで重要になるのが検出確率Cijの存在である．実際のSPECTの順問題は，図 29に示すように，
投影の他に減弱や散乱，コリメータ特性などの問題が入ってくる．その逆問題は難解となるので，でき れば順問題を確実に定式化し，順問題のみで解ける逐次近似法を用いることが望まれる．SPECTの計測 過程は図30に示すように，検出確率Cijを用いることでうまく定式化することが可能である．Cijを求め れば(22)式の繰り返しの式によって解である画像を求めることができる．実際にはCijにどこまで（減弱，
散乱，コリメータ特性）の問題を含めるかや Cijの求め方などは，メーカーや装置によって異なるが，
逐次近似法の考え方は変わらない．
最近，X線CTにも線量を少なくするなどの目的で，ML-EM法を画像再構成に用いる例が見られるよ うになった．X線CTなどの透過型CTのML-EM法は
10
1
0 ) ( 1
0 ) ( 1
0
1
0 ) ( )
( )
( ) 1 (
exp
) exp
(
N
i
M
j
k j ij i
M
j
k j ij ij N
i
i M
j
k j ij i
ij k
j k
j k j
l d
l l
y l
d l
(23)と表され，図 31 にも示すように，加減型の逐次近似法となる．加減型においてはフィードバックの値 を調整しなければならなく，(23)式においてもωなどの調整値が存在する．よって，画像を収束させる ためには，ωなどの調整値をうまく設定する必要がある．
以上のように，繰り返しの方法がどのようなことを行っているのかを，その意味合いから知ることは 重要である．実際の装置では，逐次近似法の基本式に加え，様々な工夫が盛り込まれている．その工夫 の部分を紐解くことはできないが，逐次近似法の基本をつかんでいただければ，少しでもブラックボッ クスが解消されることでしょう．特に加減型と乗除型の性質の違いだけでも感じてもらえれば幸いです．
【図の説明】
図1 順問題と逆問題の関係 図2 加減型のフィードバックの考え方
図3 加減型の逐次近似法の模式図 図4 「差分に1/20を掛けてから加える」という
フィードバックを行って繰り返した結果
図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 加算型の逐次近似法で数値シミュレ 図20 加算型の逐次近似法でα=1/900にして数値
ーションを行なった結果（α=1/1000） シミュレーションを行なった結果
図21 X線CTの乗除型の逐次近似法を 図22 乗除型の逐次近似法で数値シミュ フィードバックで示した図 レーションを行なった結果
図23 加減型と乗除型の繰り返し回数と 図24 検出確率 Cij の考え方
再構成画像の評価値
図25 16×16画素の画像と16×16投影の 図26 検出確率 Cij と投影の計算 投影データから作成した検出確率Cij
の画像（256×256画素）
図27 検出確率 Cij と逆投影の計算 図28 ML-EM法をフィードバックで表した図
図29 SPECTの順問題と逆問題 図30 検出確率 Cij による定式化
図31 透過型CTのML-EM法をフィードバック
で表した図（加減型となる）
