情報論的学習理論テクニカルレポート 2009

Technical Report on Information-Based Induction Sciences 2009 (IBIS2009)

Radon 変換を介した医用画像再構成における画像修復
Image restoration for the Medical Images using Radon Transform庄野 逸∗岡田 真人†

Hayaru Shouno

Masato Okada

Abstract: We propose an image reconstruction algorithm using Bayes inference for theRadon transformed observation data, which is usually used in medical imaging such like
CT/PET. Thorough our Bayesian reconstruction method, we introduced several hyperparameters for the prior and the observation process. The quality of the reconstructed
image is influenced of the estimation accuracy of those hyper-parameters. Hence, we also
propose the inference method of those hyper-parameters using marginal likelihood maximization principle. We show a better reconstruction result rather than that of a conventional method.
Keywords: Radon Transform, Bayes Inference, Image Reconstruction, Hyper-parameter
Inference
な画像再構成問題に対して Bayes 推定を用いた確率的画

1

まえがき
医療などで用いられる断層画像は，複数方向からの投影を行ない, 得られたデータから元の空間での信号像修復の考え方を適用することで得られる結果に関して言及する．

強度を推定することによって得られる物体断面画像で2定式化ある．このような物体のある断面における物理量の分

2.1Radon 変換

布を画像化して表示する方法を画像再構成問題 (ImageReconstruction from Projection) と呼ぶ. この画像再構成問題を計算機を通して解く場合，これは ComputedTomography (CT) と呼ばれ，現在では物体の内部断面を画像として得るために良く用いられる．CT 画像の例
としては X 線透過信号を用いた X 線 CT 画像や, 体内に注入された放射性同位体を放射線源としてその濃度分布を推定する PET (Positron Emission Tomography) 画像などが挙げられる．これらの画像再構成の問題は基本的に Radon 変換と呼ばれる投影変換を元に構成されている．本研究では，Radon 変換を観測系としてもつよう∗ 電気通信大学 電気通信学部 情報通信工学科, 〒182-8585 調布市

調布ヶ丘 1-5-1, tel. 042-443-5787, e-mail shouno@ice.uec.ac.jp,
Dept. of Information and Communication Engineering, University of Electro-Communications, Chofu-ga-oka 1-5-1, Chofu,
182-8585, Japan.
† 東京大学 大学院新領域創成科学研究科 複雑理工学専攻, 〒2778561 柏市 柏の葉 5-1-5, e-mail okada@k.u-tokyo.ac.jp,
Graduate School of Frontier Science, The University of Tokyo,
Kashiwa-no-ha 5-1-5, Kashiwa, 277-8561, Japan.

一般に Radon 変換は対象物の外部から投射される X線や，体内に注入された放射性同位体から発生する放射線等による物体を透過する信号を様々な方向から観測するような状況を想定する．Radon 変換はこのような透過してくる観測信号の強度をモデル化したものと捉えることができる．Radon 変換の概要を図 1 に示す．図 1 において観測対象の置かれた座標系を (x, y) で表すものとし，観測系の座標を (s, t) 座標で表すものとする．観測器は s 軸上にセンサアレイとして据えつけられており，観測対象を透過して落ちてくる信号を拾うことで観測投影データを構成する．このセンサアレイを観測対象の周りで 1 周回し観測を行うことで様々な方向からの投影データが得られることになる．図 1 のような状況では，センサアレイを角度 θ だけ傾けて観測しているので，(s, t) 座標系と (x, y) 座標系の間には回転変換が成立している:(st)(=cos θ− sin θsin θcos θ)(xy).

な角度 θ からの投影像に関するセンサの信号 {τ (s, θ)}yを観測し，周波数空間の任意の点の信号成分を得ること

sObject: ξ(x, y)˜ ỹ) が得られる．これをができれば，Fourier 変換像 ξ(x̃,投影定理と呼ぶ．

t

FBP 法は，投影定理を元に Fourier 逆変換 (式 (3)) を

Projection: τ(s, θ)x̃, ỹ 座標系から極座標系 s̃, θ 上での積分として変数変換θを行うことによって導出する:

x∫∫dx̃dỹ σ̃(x̃, ỹ) e2πj(xx̃+yỹ)(6)∫ π ∫ ∞=dθds̃ |s̃| σ̃(s̃ cos θ, s̃ sin θ) e2πjss̃

σ(x, y) =−∞0(7)

図 1 Radon 変換の概略: 投影データは s 軸上で観測される．θ はセンサアレイの傾きを表す．

ここで，観測対象の分布を ξ(x, y) として定義する．この分布は X 線 CT 画像の場合は位置 (x, y) における X

また式 (7) 中の s に関する逆 Fourier 変換の部分を以下のように定義しておく:∫線の吸収係数を意味し，PET 画像を再構成する場合はPositron の密度を表すことになる．観測データは，s 軸上で定義されるので，これを τ (s, θ) で表すものとすると:∫τ (s, θ) =換と呼ぶ．

2.2∫ds̃ |s̃| τ̃ (s̃, θ) e2πjss̃

g(s, θ) =dt ξ (x(s, t), y(s, t)) (1)と書けることになる．この積分変換のことを Radon 変(8)ここで∫dt ξ(x, y) =ds̃ |s̃| σ̃(s̃ cos θ, s̃ sin θ) e2πjss̃

g(s, θ) =(9)は，観測器が傾き θ のときの観測データ τ (s, θ) に対して |s̃| という微分型のフィルタをかけた像を Fouier 表現していることに他ならない．したがって |s̃| を逆 Fourier

従来手法変換した像を核関数とした畳み込み演算を行うことで，観測データ τ (s, θ) から原画像を再構成する従来手法を説明する．従来手法として良く用いられる手法として挙げられるものは Fourier 変換法および FBP (Filteredフィルタをかけた像 (Filtered Image) g(s, θ) が求められる:∫g(s, θ) =du h(u) τ (s + u),(10)Back Projection) 法と呼ばれる手法である．これらの手法はいずれも本質的には周波数空間における議論を元に

ただし, h(u) はフィルタ関数 |s̃| の逆 Fourier 変換を表理論が構築されている．本節では，FBP 法についての説すものとする．画像を再構成するには，フィルタをかけ明を行う．ここで，画像を周波数空間で扱うために，以た像 g(s, θ) に対して s = x cos θ + y sin θ の関係から下のような修復画像 σ(x, y) の Fourier 変換対を考える:∫∫σ̃(x̃, ỹ) =−2πj(xx̃+y ỹ)

dxdy σ(x, y) e(2)∫∫σ(x, y) =2πj(xx̃+y ỹ)dx̃dỹ σ̃(x̃, ỹ) e(3)一方，Radon 変換 τ (s, θ) の s に関する Fourier 変換
は, (x, y) 座標系から (s, t) 座標系への変数変換を行うことによって以下のように導ける:∫τ̃ (s̃, θ) =∫ πσ(x, y) =dθ g(x cos θ + y sin θ, θ)(11)0という積分を実行することで σ(x, y) を求めることが出来る．このように，投影データに対してフィルタリング処理を行って得た g(s, θ) から再構成を行う手法を

FBP(Filtered Back Projection) 法と呼ぶ．現在の CT
画像の再構成においても，このフィルタの畳み込み関数 h(u) をどのように設計するかということが重要な技−2πjss̃ds τ (s, θ)e(4)˜ cos θ, s̃ sin θ).
= ξ(s̃(5)術であるといえる．Ramachandran-Lakshminarayanan
らは |s̃| を逆フーリエ変換したフィルタを提案し，Logan-Shepp らは，さらに改良を加え高周波領域の抑えたよう˜ cos θ, s̃ sin θ) は，原画像の周波数空間における信号をξ(s̃なフィルタを提案している [1][2]．図 2 にこれらのフィ極座標系で表したものと考えられる．このことから，様々ルタの概形を示す．

Logan-Shepp
Frequency Domain

Ramachandran-Lakshminarayanan
Frequency Domain

●
●

●
●
●
●p(τ | σ) ∝ exp(−Hn (τ | σ))

●
●

●
●

●

●

●

●

●

●
●

●

●

●
●

●
●

●

●
●
●

●

●
●
●

●
●

●
●

●
●

●
●

●

●
●

●

-50と表す．さらに周波数空間での極座標表現を導入すると，●●●●-10(15)

●
●

●
●

●

●

●
●

●

●
●
●

●

-15とみなすことができる [3]．このとき，観測過程の確率を

●
●

●
●

510

Frequency Index

15-15-10-50510

Frequency Index

Spatial Domain

15

Spatial Domain●

Planchrel の定理より()∫∫
2p(τ | σ) ∝ exp −4π 2 γ dθ ds̃ |τ̃s̃,θ − σ̃s̃,θ |(16)

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
●

●

-15-10-505

10Posi on Index

15-15-10-5

●

●

●

●

●

●

●

●

●

●

●

●

●と書ける．ただし τ̃s̃,θ = τ̃ (s̃, θ), σ̃s̃,θ = σ̃(s̃ cos θ, s̃ sin θ)●05Posi on Index

1015である．以下では簡潔に表すために，この記法を用いることとする．

図2

古典的再構成フィルタ: 左列は Ramachandra らのフィルタ，右列は Logan Shepp らのフィルタ．上段は周波数空間での表現で，下段は実空間での表現 [1][2]．

次に Bayes 推定を行うために本研究においても画像修復を行う為の事前分布を導入する．事前分布のエネルギー関数 Hpri (σ) は観測過程の部分と同様に連続系の2.3確率モデル

積分

ここで，投影定理を元にしたモデルに対し，観測ノイズを仮定した確率モデルを定式化する．医療画像の一般∫∫dx dy ||∇σ(x, y)||2∫∫+ 4π 2 hdx dy σ(x, y)2Hpri (σ) = βのモデルとして扱う場合，放射性崩壊などのイベントは(17)

Poisson ノイズとして扱うが，本研究では可解モデルとして取り扱うことを主眼に起き，ガウスノイズを考えるのように考えるものとする．式 (17) において，第 1 項がものとしている．ここでは Radon 変換を行なう軸線上近隣の画素信号の類似度に関する拘束条件であり MRFにおける原信号 ξ(x, y) に加えてガウスノイズが重畳さに対応する．第 2 項目が画像全体の期待値の周りにどのれるものとして取り扱う．くらいの頻度で画素信号が存在するかという拘束条件で，画像上の各点においてはガウスノイズ np (x, y) を考慮

画素値全体の分布がガウス分布に従うようにするためのするので，観測値は τ (s, θ) は Radon 変換を介すること拘束条件となる．式 (17) の積分に対して，極座標表現をによって，導入したときの事前分布は∫τ (s, θ) =dt (σ(x, y) + np (x, y))(12)dtσ(x, y) + Np (s, θ)(13)∫=と 書 け る ．こ の と き Np (s, θ) は ガ ウ ス 変 数 の 積 分

Np (s, θ) =∫dt nP (x, y) となるので，やはりガウス変数として取り扱う．以上の議論と従来の画像修復のアナロ
ジーから連続系の観測過程のエネルギー関数 Hn (τ | σ)を以下のように定義する [3][5]:()2∫ π ∫∫Hn (τ | σ) = 4π 2 γdθ ds τ (s, θ) − dt σ(x, y) .
0(14)ハイパーパラメータ γ はガウスノイズ Np (s, θ) のばらつきに逆比例する量を表わし，値が大きいほどガウスノイズの影響が小さいことを意味する．これは画像修復にp(σ) ∝ exp(−Hpri (σ))(18)()∫∫= exp −4π 2 dθ ds̃ (βs̃2 + h)|s̃||σ̃s̃,θ |2(19)と書ける．
以上の議論とベイズの定理p(τ | σ)p(σ)
p(σ | τ ) = ∑σ p(τ | σ)p(σ)(20)

より事後分布を導くこととする．式 (19) および (16) を用いて式 (20) の分子の項を計算すると，()∫ π ∫2γdθ ds̃ F̃s̃ σ̃s̃,θ −p(σ | τ ) ∝ exp −4π 2τ̃s̃,θ
F̃s̃
0(())
∫ π ∫
γ2× exp −4π 2dθ ds̃ γ 1 −|τ̃s̃,θ | ,
F̃s̃
0(21)

おいて観測をガウスノイズチャネルを通して行ったものとなる. ただし，F̃s̃ = (βs̃2 + h)|s̃| + γ である．

一方，分母の分配関数∑σ p(τ | σ)p(σ) であるが，周~y波数空間を極座標に沿って離散化した上で考えるものとΔs~する (図 3 参照)．動径方向と角度方向のサンプリングΔθ(s~~k ,θl)

幅を各々 ∆s̃ , ∆θ とすると，離散化座 (s̃k , θl ) は，各々~xs̃k = k∆s̃ , θl = l∆θl と表せる．角度方向のサンプリン
グ幅は，センサアレイの観測角度に対応する．ここでは∆θ = π/Nθ とし，[0, π] を Nθ 等分した各点において観測することを仮定した．一方，動径方向のサンプリング
幅 ∆s̃ であるが，s̃ が， s の周波数表現であることから，

図3

Fourier 変換 (4) から決定する必要がある．今，長さ L散化を行う．のセンサアレイ上に Ns 個のセンサを等間隔に配置することを考える．このとき，Fourier 変換 (4) を離散化して表現すると，∫τ̃ (s̃k , θl ) =∼周波数空間における離散化方法: 極座標上での離考慮すると

∫ π

∫ ∞

hσ(x, y)i =dθ−∞

0ds τ (s, θ)e−2πjss̃N∑s −1(22)∆s τ (sk , θl )e−2πjsk s̃k̃(23)k̃=0と書ける．このときのナイキスト周波数は，1/2∆s であり，対応する離散フーリエ変換は，周波数空間の定義域

ら，周波数空間のサンプリング間隔は，∆s̃ = 1/Ns ∆s とFourier 空間における空間周波数の集合 {hσ̃s̃,θ i} を求めれば良く， hσ̃s̃,θ i は形式的に∑hσ̃s̃,θ i =σ σ̃s̃,θτk,l e−2πj Ns = ∆s τ̃k̃,lkk̃

(24)k̃=0skk̃
1 ∑τ̃k̃,l e−2πj NsNs(25)k=0という離散フーリエ変換対 {{τk,l }, {τ̃k̃,l }} を用いて表現できる．ここで，式 (20) に含まれる 積分∑Ns −1k̃∫ds̃ を∆s̃ として離散化を行うと，周辺化事後確率)(Nsγ
τ̃k̃,l , 2(26)p(σk̃,l | τ ) = N σk̃,l |Fk̃
8π ∆θ Fk̃が得られる．ただし，N (x | µ, S) は，x を確率変数とした平均 µ, 分散 S の正規分布の確率密度を表すものとする．また，周波数空間上の各サンプリングポイント

σ̃k̃,l τ̃k̃,l , F̃k̃ は，それぞれ σ̃k̃,l = σ̃s̃k̃ ,θl , τ̃k̃,l = τ̃s̃k̃ ,θl ,F̃k̃ = F̃s̃k̃ = (βs̃2k̃ + h)|s̃k̃ | + γ を表すものとした．

2.4

画像修復

画像修復の戦略としては，周辺事後分布による期待値hσ(x, y)i を取る戦略を考える．Fourier 表現の線形性を

0dθ

R∞−∞

˛
˛2˛
˛ds̃ F̃s̃ ˛σ̃s̃,θ − F̃γ τ̃s̃,θ ˛s̃˛˛2˛˛dθ −∞ ds̃ F̃s̃ ˛σ̃s̃,θ − F̃γ τ̃s̃,θ ˛0

RπR∞s̃

(28)と書ける．式 (28) を計算する為に前節と同様に周波数空間の極座標上で離散化を行なって解を求めると

DEγτ̃σ̃k̃,l =F̃k̃ k̃,l(29)となる．これは式 (26) において MAP (Maximum APriori) 戦略を採った時の解と一致する．
{DE}
以上により求まったN −1

τ (sk , θl ) ∼e−4π 2σe換対は，
N∑s −1Rπ−4π 2∑

書けることになる．以上の議論から式 (4) の Fourier 変τ̃ (s̃k , θl ) ∼ ∆s(27)と 書 け る ．し た が っ て ，hσ(x, y)i を 求 め る た め に は,を [−1/2∆s , 1/2∆s ] とし，この区間を Ns 等分したサンプリング点での値で表現されることになる．このことかds̃ |s̃| hσ̃s̃,θ i e2πjs̃(x cos θ+y sin θ)σ̃k̃,lを用いて，FBP 法に出てくる関数 g(s, θ) の位置 (sk , θl ) における期待値

hgk,l i = hg(sk , θl )i を計算すると∫ ∞hgk,l i =ds̃ |s̃| hσ̃s̃,θl i e2πjsk s̃−∞sDEkk̃
1 ∑|s̃k | σ̃k̃,l e2πj Ns
Ns

N −1

∼

(30)

k̃=0となる．
この hgk,l i を元に画像空間での位置 (x, y) に対応する画素値の平均 hσ(x, y)i を推定することを考える．この場合は s = x cos θ + y sin θ という関係があるので，s が丁度サンプル格子点に乗るとは限らない．そこで線形補間で画素値を推定することにする．すなわち sk∗ ≤ s <k∗としたsk∗ +1 となるような，k ∗ を探し出し，v = s−s∆s上で，hg(s, θl )i = (1 − v) hgk∗ ,l i + v hgk∗ +1,l i(31)表1

Histogram of pixels

3.0 x 104

推定したハイパーパラメータ値．事前分布に関わるハイパーパラメータ β, h は画像の事前分布 (式 (17)参照) に，γ は観測過程 (式 (14) 参照) に関わるハイ

2.0 x 104

パーパラメータ．
1.0 x 104

標準偏差

βh−3

γ−2

0.5

8.54 × 10

3.28 × 10

1.65 × 104

1.0

7.80 × 10−3

3.62 × 10−2

2.76 × 103

図 4 Logan-Shepp ファントム画像と画素値のヒスト

2.07.14 × 10−3

3.94 × 10−2

4.96 × 102

グラム．Logan-Shepp ファントム画像は体断面を模式

4.06.73 × 10−34.17 × 10−2

1.07 × 102

0123456的に表現し，CT/PET 画像の再構成シミュレーションとして用いられる．図はプログラミング言語 R の拡張表現できる:パッケージによって生成している．右図はこのファン

N −1

トム画像の画素値のヒストグラムを表す．背景色の黒ln Zpri (β, h) = −は画素値 0 で表され，もっとも明るい白色も 6.0 程度の画素値で表現されている．と 置 く こ と に す る ．最 終 的 に 求 め た い 再 構 成 の 位 置(x, y) における期待値 hσ(x, y)i は∫ πhσ(x, y)i =dθ hg(x cos θ + y sin θ, θ)i∼ ∆θNθ Nsln γ(39)2)
Ns −1 (24π 2 ∆θ ∑γ
ln Zpost (β, h, γ) = −γ 1−τk̃,lNs
F̃kln Zn (γ) = −k̃=0(32)−hg(x cos θl + y sin θl , θl )i(33)Nθ2N∑
s −1ln F̃k̃ .

(40)

k̃=0

l=0

あとは，周辺対数尤度を各ハイパーパラメータに関しと近似できる．

2.5(38)k̃=00N∑θ −1sNθ ∑ln(βs̃2k̃ + h)
2て勾配を求め，これを 0 とおくことによってハイパーパラメータ推定を行うことができる．ただし解析的に解く

ハイパーパラメータ推定ことは困難であるため，ハイパーパラメータ推定は，式画像修復を行なうためには，ハイパーパラメータ β, h,γ を推定する必要がある．これらのハイパーパラメータ(34) を最大化するような，β, h, γ を勾配法で求めることとした．すなわちを対数周辺尤度最大化の枠組みで求めることを考える．

  t β t+1β ht+1  =  ht  + η∇ ln p(τ | β, h, γ)
γ t+1γt

対数周辺尤度はln p(τ | β, h, γ) = ln Zpost (β, h, γ)− ln Zn (γ) − ln Zpri (β, h).

(34)のように 3 つの分配関数，すなわち事後分布の分配関数

Zpost , 観測チャネルの分配関数 Zn , 事前分布の分配関数
Zpri ，各々の対数で記述でき，これらの分配関数は各々∑Zpri (β, h) =e−Hpri (σ | β,h)(35)
σZn (γ) =∑e−Hn (τ | σ,γ)(36)

e−Hpri (σ | β,h)−Hn (τ | σ,γ)(37)

τZpost (β, h, γ) =∑σ(41)によってハイパーパラメータを更新する．

3計算機実験
計算機実験における座標系は以下のような形で取るものとした．まず，Logan-Shepp ファントム画像を Nx ×Ny の大きさで生成し，これを原点を中心とした縦横 Lの長さの正方形 [−L/2, L/2] × [−L/2, L/2] 内の領域にマッピングしたものを原画像とした．また原点を中心とした半径 L/2 の円の外側はセンサアレイに投影したときに, 観測角度によっては観測できない場合もあるため，とかける．各分配関数の積分に対して前節までの議論とこの部分には観測対象は存在しないものとした．また各同様な離散化を行うと，定数部分を除いて以下のように観測角度 θl において s の座標系は，図 1 のように考え，原点 s = 0 の位置は，(x, y) 座標系の原点を s 軸上に垂直に投影した点にとるものとした．今回の実験で

FBP
Reconstrucon

Bayes
Reconstrucon

Noise S.D. = 0
図6Noise S.D. = 1.0

Noise S.D. = 2.0

Noise S.D. = 4.0

再構成画像の画質比較． 上段は FBP 法による再構成結果，下段は Bayes 法でハイパーパラメータ推定を行った後に再構成した結果. 各列は観測時に加えるノイズの標準偏差を変えて再構成した結果を示している．

PSNR [dB]
40
●

Bayes
FBP
●

●

30●
●
●
●
●

201.0

図72.0
3.0
Noise S.D.

4.0

再構成画像の画質の数値比較．横軸は観測時に加えるノイズの標準偏差を表し，縦軸は PSNR を示している．PSNR 値は高いほうが画質が良いことを表す．

は， Nx = Ny = Nθ = Ns = 256 とし，センサの長さは

L = 1 とした. 図 4 に今回の実験で用いた Logan-Shepp
図5推定された再構成フィルタの周波数空間での表現．加えるノイズの標準偏差の大きさを表 1 に従って設定した．標準偏差の大きさは，上段左からは 0.5, 1.0,下段左からは 2.0, 4.0 に各々設定した．ノイズの大きファントム画像と，その画素値のヒストグラムを示す．
これからわかるように原画像の画素値は [0, 6] の値で表現されている．

さが小さい時は図 2 左の様なフィルタを構成するのに観測過程においては，Radon 変換においてガウスノイ対して，ノイズが大きくなると，図 2 右の様なフィルタ

ズ np (x, y) が加わるものとした．このとき np (x, y) の標を f 構成し，さらにノイズが大きくなると大きく高周波成分が抑制される．

準偏差が，0 〜 4 の範囲で変化させ，ハイパーパラメー
タ推定を行った画像の画質評価を行った．
表 1 に式 (41) によって推定したハイパーパラメータ値を示す．画像の事前分布に関わるハイパーパラメータ

β, および h の大きさはほとんど変化していないのに対離散化を極座標系上で行うことによって，事後分布が可して，観測過程に関わるハイパーパラメータは，ノイズ解モデルとして取り扱うことができるようになることをの標準偏差が大きくなるほど観測精度を示す γ の値が小示し，これを用いて MAP 推定値および，平均ピクセルさくなっていることから，定性的には妥当な結果である

値 hσ(x, y)i の推定の手続きを記述した．

ように見える．また，それぞれの推定ハイパーパラメー

さらに事後分布に含まれるハイパーパラメータ β, h, γ

タ結果に対して，得られた再構成フィルタの周波数空間を周辺対数尤度最大化原理から求める手続きを示し，計での表現を図 5 に示す．図よりノイズの標準偏差が小さ算機実験としてハイパーパラメータを求めた上で画像再い時は，Ramachandran-Lakshminarayanan 型のフィル構成を行った結果を示した．

タ関数 (図 2 左上参照) であるが，重畳されるノイズの標準偏差が大きくなるにつれ，高周波成分を落とすような参考文献フィルタが自動的に構成されていることがわかる．

[1] G. N. Ramachandran and A. V. Lakshmi-

これらの推定ハイパーパラメータ値を用いて画像再narayanan: Three-dimensional reconstruction from

構成の計算機シミュレーションを行った．得られた再radiographs and electron micrographs, Proc. Natl.

構成結果を 図 6 に示す．図 6 の上段は従来手法である

Acad. Sci., vol.68, pp.2236–2240, 1971.

FBP 法 を用いて再構成を行った結果である．比較対象[2] L. A. Shepp and B. F. Shepp: The Fourier recon-の FBP 法としては数理的なアルゴリズムとして導かれ

struction of a head section, IEEE Trans. Nucl. Sci.,

る Ramachandran-Lakshminarayanan らの手法を用い

vol.NS-21, pp.21-43, 1974.ている [1]．一方，図 6 下段は本研究で提案しているハイ[3] K. Tanaka and J. Inoue: Maximum likelihood hy-

パーパラメータ推定を行った上で Bayes 手法を用いて再perprameter estimation for solvable Markov ran-構成を実行した結果を示している．目視においても従来

dom field model in image restoration, IEICE手法である FBP 法においては混入するノイズが大きく

Trams. Inf. Syst., vol. E85-D, No.3, 2002.

なるほど，画質が劣化していくのが目視で確認できる．

[4] K. Tanaka: Statistical-mechanical approach to im-

さらに画質の数値比較を行う手段として PSNR (Peakage processing, J. Phys. A, vol. 35, No.37, pp. 546

Singal to Noise Ratio) を用いて画質評価を行った結果– 557, 2002.

を，図 7 に示す．図 7 の横軸は観測時に加えるノイズの標[5] K. Tanaka, H. Shouno, M. Okada and D. M. Titter-

準偏差を表し，縦軸には PSNR 値を表している．PSNR

ington: Accuracy of the Bethe Approximation for

値は値が大きいほど画質劣化が少ないことを表す．図中

Hyperparameter Estimation in Probabilistic Imageのグラフからわかるように Bayes 手法を用いた場合は，Processing, J. Phys. A, vol. 37, No.36, pp. 8675 –PSNR の意味において再構成時における画質劣化が抑え

8696, 2004. pp.329-338, 2003.

られ，ノイズの標準偏差が 4.0 の場合にでも 27.5 [dB] 程度を維持できることが示されている．一方 FBP 法では標準偏差が 1.5 程度で，既に PSNR 値が 27.7 [dB] 程度まで下がってしまう．以上の議論から，周波数フィルタに工夫をしていない FBP 法と比べて，我々の提案する手法の方が観測ノイズに関して耐性が大きいということがいえる．

4

まとめRadon 変換を用いるような観測系における画像修復の議論を Bayes 推定の枠組みで行なった．事前分布として
は，隣接ピクセル間での滑らかさを制御する MRF 的なガウス分布 p(σ) を導入し，観測過程はガウスチャネルを通して行なわれるものとして，p(τ | σ) を定式化した．