研究要旨 - 室蘭工業大学

強度分布特性を活かした進化型多目的最適化に基づく
少数投影 CT 再構成アプローチの検討
長舟 和馬 † , 渡邉 真也 †† , 塩谷 浩之 ††
‡
室蘭工業大学 情報電子工学系学科 † , 室蘭工業大学大学院 しくみ情報系領域 ††
1
はじめに
一 般 的 な コ ン ピュー タ 断 層 撮 影 (Computed
Tomography, CT) では対象となる物体 (原画像) に
対して 360 度の全方向から X 線を透過し, その分
布を利用することで断面画像を再構成してい
る 1) .し か し, 障 害 物 な ど に よ り 投 影 方 向 に 制
限が存在し少数でしか投影が行えない状況もあ
り, 一 般 に 少 数 投 影 CT(Sparse CT) と 呼 ば れ る .
この問題は不完全問題の 1 種であり, 原画像を復
元するためには少数方向からの観測情報, 物体の
存在位置といった一部の情報を手がかりに断面
画 像 を 推 定す る必 要 がある 2) .
現実問題として, 欠損情報が多量なればなるほ
ど補完に対する計算負荷は増大するが, 近年の計
算機環境の飛躍的向上に伴い大いに盛り上がり
を 見 せ て いる 3, 4) .
少 数 投 影 CT に 対 す る 先 行 手 法 と し て, フィル
タ補正逆投影法 (Filtered Back Projection 法, FBP
法) や空間適応フィルタの反復適用による再構成
(Recursive Spatially Adaptive Filtering, RSAF) な
ど複数の手法が提案されているが, それぞれ幾つ
か の 問 題 点が 指 摘 され て いる 5, 6) .
そこで, 本研究では先行研究 7) である, 少数投影
CTに対して進化型多目的最適化手法(Evolutionary
Multi-criterion Optimization, EMO) 8) の 枠 組 み
を利用した画像再構成手法に対して強度分布特
性を活かした遺伝的操作を新たに考案し, 従来研
究 に 対 す る推 定 画 像の 精 度向 上 を 試 み た .
An evolutionary multi-criterion optimization approach utilizing the characteristics of strength distribution for sparse CT image reconstruction
†
Kazuma NAGAFUNE(s2124122@mmm.muroranit.ac.jp)
††
Shinya WATANABE(sin@csse.muroran-it.ac.jp)
‡
Hiroyuki SHIOYA(shioya@csse.muroran-it.ac.jp)
Department of Infomation and Electronic Engineering,Muroran Institute of Technology (†)
College of Information and Systems, Muroran Institute of
Technology (††)
本問題において推定する対象は実空間, 周波数
空間の両空間であり, それらを同時に最適化する
必 要 が あ る た め, こ こ で は 多 目 的 最 適 化 問 題 と
し て 定 式 化 し た .さ ら に, 本 問 題 に お い て 推 定
すべき欠損情報は画素数と比例し膨大である
ため, 本手法では大域の探索では EMO を使い多
点探索を行ない, 局所探索では GS アルゴリズム
(Gerchberg-Saxton algorighm) 9, 10) を利用したハ
イブリッド方式を採用した.
また, 本論文では遺伝的操作として, 周波数空
間における上下左右の対称関係に着目した新た
な交叉, 突然変異手法について提案, 有効性の検
証を行った.具体的には, 周波数空間の強度分布
では 4 分の 3 に強く類似性があることが知られて
いるため, この特性を陽に遺伝的操作に取り入れ
ることによる探索効率の向上を試みた.
本提案手法は, 多点に基づく多目的最適化探索
であるため一度の探索で特徴の異なる多様な解
候補を導出することが可能である.これは, 本対
象問題が逆問題の 1 種であり一意に解を求めるこ
とができないという特性に対して非常に効果的
であると考えている.
以下, 少数投影 CT についての概説, 提案アルゴ
リズムの流れと強度分布特性を生かした交叉, 突
然変異手法を紹介し, 数値実験と結果, 結論を述
べ る.
CT の原理と少数投影 CT
2
本節では, 少数投影 CT による画像再構成につ
い ての概要と先行手法の紹介を行う.
2.1
CT の基本原理
コンピュータ断層撮影 (ComputedTomography,
CT) は全方向からの投影による透過量を計測す
ることで断面画像 † を得る撮影方法であり, 複数
の 再 構 成 手 法 が 提 案 さ れ て い る 1) .代 表 的 な
†
対象物からの透過量の吸収分布を投影データ(Projection
Data) と 呼 ぶ .
3.2
Fig. 1 The relationship Object and Inverse domains based on Fourier transform.
例 と し て, 計 算 機 を 必 要 と し な い 単 純 重 ね あ わ
せ 法, フ ー リ エ 変 換 を 利 用 す る 解 析 的 手 法 で あ
る フィル タ 逆 投 影 法 (Filtered Back-Projection 法,
FBP 法) やそれに類似した重畳積分法など複数の
手 法 な ど が挙 げら れ る 6) .
2.2
少数 投 影 CT
少数投影 CT は, ごく限られた方向からの投影
データを基に画像再構成を試みる問題であり, 不
完全問題の 1 種である.少数投影 CT 画像では, 投
影方向数が少数であればあるほど画像の解像度
に比例して欠損情報量も増加するため, 原画像の
推 定 が 指 数的 に 困 難に な ると い う 特 徴 を 持つ.
少数投影 CT における先行手法として RSAF と
呼ばれる空間適応フィルタの反復適用による再
構成があり, 2.1 節で紹介した FBP 法と合わせて
次 章 で 解 説す る.
少数投影 CT に対する先行研究
3
本節では 2.1 節と 2.2 節で紹介した先行手法の内
FBP 法, RSAF に つ い て解説 す る .
3.1
フィルター 逆 投 影法
フィル タ 逆 投 影 法 (Filtered Back-Projection 法,
FBP 法) は, Fig. 1(a), (b) にあるように, 指定した
方向からの実空間の投影データを周波数空間に
対 応 す る 角 度 に 挿 入 し, 得 ら れ た 周 波 数 空 間 の
データをフィルタリングの後に逆投影すること
で 画 像 再 構成 を 行 なって いる.
ここで, (x, y) 座標系をある角度 θ 回転させた座
標系を (r, s) とするとき, f (x, y) の θ の投影データ
h(r, θ) は以 下の 式 (1) で与 えられ る .
∫
f (r cos θ − s sin θ, r sin θ + s cos θ)ds (1)
h(r, θ) =
R
空間適応フィルタの反復適応による再構成
観測対象の信号のスパース性や圧縮可能性の
仮定のもとでその少数データだけから観測信号
を復元するための技術は, Compressed Sensing(も
しくは Compressive Sampling) と呼ばれ, 広く研究
が行われている 5) .
本研究では空間適応フィルターの反復適用によ
り推定画像の修復と更新を実現する Egiazarian ら
の提案する手法 11) (Recursive Spatially Adaptive
Filtering, RSAF) を 採 り 上 げ る .RSAF は 実 空 間
上での空間適応フィルター Φ を利用することで観
測データのスパース性を補完する.Φ の適用によ
る画像更新に基づいて周波数空間未知領域の
推 定, 実 空 間 画 像 の 導 出 を 行って い る .な お, Φ
として Block-Matching and 3D filtering algorithm
12) を使用している.
EMO に 基 づ く 少 数 投 影 CT 再 構 成 へ の
アプローチ
4
少 数 投 影 CT に 対 し, EMO と GS ア ル ゴ リ ズ ム
(Gerchberg-Saxton algorithm)9, 10) を組み合わせ
た再構成手法の利点と全体の流れを述べ る.
提案手法の利点
EMO による多点探索という特徴と GS アルゴリ
ズムにおける拘束条件を利用した効率の良い探索
を組み合わせることによる主な利点を下記に示す.
• 多点探索により一度の試行で特徴の異なる
複数の解 候補を導出可能
• ブ ロック マッチ ン グ 法 の よ う に 実 画 像 に 関
する情報 を全く利用せずに推定可能
提案アルゴリズムの流れ
Fig. 2 に示す本提案手法の具体的な手順を以下
に 示す.
Step 0 EMOの各種パラメータの初期化と初期個
体・母集団の生成と個体の適合度を評価
Step 1 EMO による探索
Step 2 各個体に対する GS アルゴリズムの適用
Step 3 終了条件判定 (条件を満たしている場合
は終了, それ以外は Step1 へ)
4.1
周波数空間の強度分布特性を活かした 突
然変異と交叉
本 研 究 で 対 象 と し て い る 少 数 投 影 CT 問 題 は,
周波数空間における未知パラメータの推定に帰
⇒⇒
⇒⇒
⇒⇒
similar
similar
similar
Fig. 4 The concept figure of mutation method.
Fig. 2 Flowchart of image reconstraction of
Sparse CT based on EMO.
⇒⇒
⇒⇒
⇒⇒
similar
similar
similar
Fig. 3 Symmetric property of intensity distribution in frequency space.
い, それらの重みづけ和の値により各領域を更新
する.この突然変異を利用することにより, 4 つ
の領域の類似性が向上し探索する設計空間を疑
似的に削減する効果が見込まれるため, 大幅な探
索効率の向上を期待することができる.
ただし, 厳密には Fig. 3 1 や Fig. 3 2 間に示さ
れるように各領域は完全に一致する対応関係に
はない.そのため, Fig. 4 の太い矢印に示すよう
に自身の領域の重みを大きく, それ以外の 3 つの
領域部分の重みを軽くすることにより, その領域
独自の特徴がある程度残るような方法を用いた.
具体的には, 自身の領域の重みを ωi , それ以外の
領域の重みを ωj {j = 1, 2, 3} として以下式 (3) に示
す更新式を突然変異として用いた.
Gi (x, y)
=
Gi (x, y)ωi +
∑
Gj (x, y)ωj (2)
j∈I\{i}
着する.そのため, 本研究では, 周波数空間にお
ける実部と虚部を実数値で表現したものを設計
パラメータ(遺伝子)として扱い, 各種遺伝的操
作 を 行って い る.
周波数空間における強度は, Fig. 3 に示すよう
に一定の対称性を有することが知られているた
め, 遺伝的操作においてこの強度分布特性を陽に
利用することにより, さらなる探索効率の向上を
期待することができる.以下, 本研究で新たに提
案 す る 突 然変 異, 交叉 手 法につ い て 説 明 す る.
突然 変 異 手 法
周波数空間における強度の対称性から 1 つの画
像は Fig. 3 に示すように 4 つの領域に分割するこ
と が で き, そ れ ぞ れ の 領 域 間 は 強 い 類 似 性 を 有
している.そこで本研究では, 突然変異の中に 4
つの領域の類似性を向上させる仕組みを陽に
組み込む手法を考案した.具体的には, 1 つの画
像 を Fig. 3 に 示 す 4 領 域 に 分 割 し, 各 領 域 に 対 し
て Gi (x, y), Gi (x, y){i = 1, 2, 3, 4} の重みづけを行
s.t.
s.t.
交叉手法
0.5 ≤ ωi < 1.0
∑
ωj = 1.0
ωj +
(3)
(4)
j∈I\{i}
本手法に求められる重要な要素の 1 つは, 特徴
の異なる多様な解を生成することである.その
ため, ここでは 2 つの親から 2 つの子を生成する典
型的な交叉ではなく, 複数の親から複数子を生成
する方法を採用した.
提 案 す る 交 叉 手 法 で は, 突 然 変 異 手 法 と 同 様,
Fig. 3 に あ る 対 称 性 に 基 づ き 4 分 の 1 の 領 域 ご と
に画像を区分けし, 各領域の類似性を向上させる
仕組みに基づいている.ただし, 複数の親の形質
を反映させるため, 突然変異手法と同様に各領域
の重みづけ和を各親ごとに求め, それらの値を合
算, 平均化した値により子のパラメータ値を算出
している.
具 体 的 に は, ま ず 式 (3) と 同 様 の 方 法 に よ り 各
親ごとに重みづけ和を求め, 次にその値について
全ての親の平均値を算出し, それを子どものパラ
メータ値とした.この操作を生成する子供の数
Separate
Paren1
Separate
Paren2
Fig. 6 Real constraint condition and picture image of LENA.
Combine
min f1
= Eout (g) + Eimage (g) + Epos (g)
(5)
min f2
= Ediff (G)
(6)
Eout (g) :実空間拘束条件の範囲における不正値
の合算値
Kid1
Fig. 5 The concept figure of crossover.
Eimage (g) :実空間拘束条件の範囲外における負
の値を絶 対値化した合算値
Table 1 Parameters of NSGA-II.
Epos (g) :実空間の虚部における, 0 以外の値の合
算値
Parameter name
Number of direction of projection
Number of population
Number of GS algorithm
Terminal generation
Number of trials
Crossover rate
Mutation rate
Value
4, 8, 16
40
20
10
10
1.00
0.004
だけ繰り返すことにより, 多くの親から多くの子
どもを生成する交叉を実現している.本交叉手
法 の 概 念 図を Fig. 5 に 示す.
5
数値実験
本手法の有効性を検証するため, 解像度256×256px
で複雑度の異なる Phantom 画像, LENA 画像 (Fig.
7, Fig. 6) の 2 種類に対する数値実験を行った.実
験 で は, 投 影 数 が 4, 8, 16 方 向 の 3 つ の 場 合(Fig.
8)に お け る, FBP 法, RSAF, 提 案 手 法 の 結 果 を
比 較 し た.
5.1
実験 設 定
提案手法における EMO アルゴリズムとしては,
実数値型 NSGA-II 13) を利用し, 交叉手法と突然
手 法 は 4.1 節 で 提 案 し た 手 法 を 用 い た .ま た, 本
手法である EMO アルゴリズムの設定パラメータ
は, 個体数を 40 個とし, 世代ごとに行なう GS アル
ゴリズムの回数を 20 回とした.その他使用パラ
メ ー タ を含 めた 一 覧 を Table 1 に 示 す.
ここでは, 下記に示す 4 種類の評価基準を各空
間ごとに合成することで 2 目的最適化問題として
定 式 化 し た.
Ediff (G) :周 波 数 空 間 に お け る 観 測 デ ー タ の 誤
差の合算 値
な お, RSAF の 繰 り 返 し 設 定 に つ い て は, 提 案
手法で利用している計算量の多いフーリエ変換
と逆フーリエ変換を行なうサイクルの数だけ設
定した.具体的には, 1 回の試行で画像出力に関
するサイクル数は 60 回, アルゴリズム内で行われ
る サ イ ク ル 数 は 400 回 で あ る た め 合 計 460 回 を
RSAF の 繰 り 返 し 数 と し た .一 方, RSAF に お い
て推奨される繰り返しは 10000 回であるため, そ
の 数値における実験も行なった.
5.2
評価方法
本研究では解像度M × N の原画像f と推定画
像 g の異なり具合を定量的に評価するため, 推定
した実画像の評価として式 (8) に示す PSNR を設
定 し た .PSNR は 一 般 的 に 使 用 さ れ る PSNR と
同じであり, 最大値は最大輝度 (今回はグレース
ケールのため 255) である.
M N
1 ∑∑
MSE =
{f (i, j) − g(i, j)}2 ,
M N i=1 j=1
)
(
MAX2
.
PSNR = 10 · log10
MSE
5.3
(7)
(8)
実験結果
各手法により得られた推定画像をそれぞれ Fig.
9, Fig. 10 に示す.なお, 参考として RSAF の推奨
パラメータによる実験結果を Fig. 11 に示 す.
Fig. 7 Real constraint condition and picture image of Phantom.
Fig. 8 Projection direction: 4, 8, 16 radial lines.
Fig. 11 Applying RSAF by recommended iterations.
Table 2 The results of PSNR values for Phantom.
Each method
FBP
RSAF
The proposed
(best results)
The proposed
(average results)
FBP
RSAF
Proposed
Fig. 9 Result of reconstruction in Phantom picture images at radial 4,8,16 lines.
4 lines
11.341
15.258
Num of projections
8 lines 16 lines
11.652
12.126
15.759
16.388
13.733
14.577
13.739
14.583
16.925
16.930
ま た, 各 手 法 に お け る PSNR 値 に つ い て Table
2, Table 3 に示す.なお, 提案手法による PSNR 値
は, 全試行の中で最も PSNR 値の高かったものと,
その試行における全個体の平均をとったものを
示し, 全個体平均がどの個体よりも良好な画像と
なっていることを確認する.
5.4
考察
Fig. 9 お よ び Fig. 10 か ら, 投 影 方 向 数 が 増 加
するほど得られる実像の精度が向上することが
分かる.また, PSNR 値の高い母集団の個体平均
は, どの個体よりも高い評価値の推定画像となっ
ていることも確認することができる.このこと
は, 少数投影 CT において一意な解は存在せず, 特
徴の異なる多様な解の平均を求めることが, 結果
Table 3 The results of PSNR values for LENA.
Each method
FBP
RSAF
Proposed
Fig. 10 Result of reconstruction in LENA picture
images at radial 4,8,16 lines.
FBP
RSAF
The proposed
(best results)
The proposed
(average results)
4 lines
09.933
15.972
Num of projections
8 lines 16 lines
10.094
10.235
18.887
21.388
13.553
15.203
15.821
13.558
15.213
15.838
として原画像に最も近い推定となっていること
を 示唆 し て い る .
ま た, 先 行 研 究 に お け る 補 完 機 能 を 持 た な い
FBP 法では, RSAF および提案手法と比較して適
切に再構成が行えていないことが分かる.この
ことより, 原画像の推定のためには欠損情報に対
する何らかの補完が必須であることが分かる.
一方 RSAF は, たとえ繰り返し数が推奨値より十
分に小さい場合でも提案手法に比べ優良な結果
が 得 ら れ て い る こ と が Table 2, Table 3 か ら 読
み 取 れ る .こ れ は, 少 数 投 影 CT に お い て ブ ロッ
ク マッチ ン グ 法 を 実 画 像 フィル タ ー と し て 用 い
ることの有用性を示しており, 実画像データの持
つ特徴に着目した欠損部分の推定が少数投影 CT
に お い て も効 果 的 であ る こと が 分 か る .
6
おわりに
本 論 文 で は 少 数 投 影 CT に 対 し, EMO と GS ア
ルゴリズムの組み合わせたアプローチの有効性
に つ い て 検 証 を 行った .提 案 手 法 の 大 き な 特 徴
は, 周波数空間の特性を遺伝的操作に陽に取り入
れている点であり, 紙面の都合上示さなかったが,
一般的な SBX や突然変異を使用した場合に比べて
探索効率を大幅に向上させた手法となっている.
本提案手法は, 実空間画像における特徴を陽に
扱 い 欠 損 情 報 の 補 完 を 行 う RSAF な ど と は 根 本
的に異なったアプローチをとっており, 計測デー
タの欠損部分を GS アルゴリズムの概念を利用し
て 愚 直 に 推定 す る もの で ある .
2 種類の複雑度の異なる画像を用いて数値実験
を行った結果, 以下の事柄を明らかにすることが
できた.
• 周波数特性を考慮した遺伝的操作を組み込
むこ とに よる 有 効性
• 対称画像が単純である場合には, 欠損部分が
多い場合においても真の解に近い画像を導
出可能
• 特徴の異なる複数の推定画像の平均化は, 真
の 画 像 に 近 い良質 な 画像と な る .
一方, Lena 画像のように複雑において欠損情報
が多い場合には, 真の解に近づくことが困難であ
る こと な ど が 明 らか と なった .
今後は, 単純比較では劣っていた RSAF の特徴
である画像フィルタリング手法について調査し,
提案手法との組み合わせについて検討を進めた
いと考えている.
参考文献
1) 日本放射線技師会. CT システム入門コンピュー
タ断層撮影の理論と実際. マグブロス出版, 1991.
2) R. Gordon. A tutorial on art(algebraic reconstruction techniques). IEEE Trans-actions on Nuclear
Science, 1974.
3) 岩間尚文. 少数投影計算機トモグラフィーの最適
設 計. 文 部 省, 1994.
4) 矢島邦昭, 田山典男. 三角格子標本化を用いた少
数投影データからの海洋音速分布画像の再構成
シミュレーション. 計測自動制御学会東北支部第
221 回 研 究 集 会, 2005.
5) 平 林 晃. Compressed sensing: 基 本 原 理 と 最 新 研
究 動 向. 電 子 情 報 通 信 学 会 技 術 研 究 報 告. VLD,
VLSI 設 計 技 術, 2009.
6) 篠原広行. エクセルによる画像再構成入門. 医療
科 学 社, 2007.
7) 大江将吾. 進化型多目的最適化に基づく少数方向
投影からの CT 画像再構成,室蘭工業大学,修士
論 文. 2011.
8) K. Deb. Multi-objective optimization using evolutionary algorithms. Wiley, 2001.
9) R. W. Gerchberg and W. O. Saxton. A practical algorithm for the determination of phase from image
and diffraction plane pictures. Optik, 1972.
10) 塩谷浩之, 郷原一寿. 位相回復–計算アルゴリズム
(ミニ特集回折イメージング–位相回復の新展開).
計 測 と 制 御, 2011.
11) K. Egiazarian, A. Foi, and V. Katkovnik. Compressed sensing image reconstruction via recursive
spatially adaptive filtering. In Image Processing,
2007. ICIP 2007. IEEE International Conference
on, 2007.
12) K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian. Image denoising with block-matching and 3D
filtering. In Proceedings of SPIE, 2006.
13) 渡邉真也. 近傍個体の交叉に基づく多目的遺伝的
アルゴ リズム とその 応用に 関する研 究,博士 論
文 ,同 志 社 大 学. 2003.