強度分布特性を活かした進化型多目的最適化に基づく少数投影 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([email protected]) †† Shinya WATANABE([email protected]) ‡ Hiroyuki SHIOYA([email protected]) 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 = Ediﬀ (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 目的最適化問題として定式化した． Ediﬀ (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 diﬀraction 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.