Comments
Description
Transcript
適応マスクを用いたトランケート劣化画像の ブラインド回復と応用 石原信人
適応マスクを用いたトランケート劣化画像の ブラインド回復と応用 Blind Image Recovery of Truncated Degraded Images Using Adaptive Mask and Its Applications 2008 年 2 月 石原 信人 適応マスクを用いたトランケート劣化画像の ブラインド回復と応用 Blind Image Recovery of Truncated Degraded Images Using Adaptive Mask and Its Applications 2008 年 2 月 早稲田大学大学院 理工学研究科 物理学及応用物理学専攻 光物理工学研究 石原 信人 iii 目次 第1章 序論 1 第2章 原理 5 2.1 画像処理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 デコンボリューションフィルター . . . . . . . . . . . . . . . . . . . . . 8 2.3 ブラインド・デコンボリューション . . . . . . . . . . . . . . . . . . . . 9 2.4 まとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 最適化法 13 最適化法とは . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 第3章 3.1 3.2 3.1.1 最急降下法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.1.2 山登り法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.1.3 共役勾配法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.1.4 最尤法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.1.5 遺伝アルゴリズム . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.1.6 ニューラルネットワーク . . . . . . . . . . . . . . . . . . . . . . 16 3.1.7 焼き鈍し法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 ブラインド・デコンボリューションの最適化法 . . . . . . . . . . . . . . 17 3.2.1 フーリエ反復法 (AD アルゴリズム) . . . . . . . . . . . . . . . . 18 3.2.2 ゼロシート法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 iv 目次 3.2.3 焼き鈍し法 (SA 法) . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2.4 ハイブリッド AD/SA 法 . . . . . . . . . . . . . . . . . . . . . . 24 3.2.5 その他のアルゴリズム . . . . . . . . . . . . . . . . . . . . . . . . 25 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.3.1 原理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.3.2 結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 まとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 適応マスク手法 31 4.1 従来手法のブラインド・デコンボリューションの問題点 . . . . . . . . . 31 4.2 トランケートが画像に及ぼす影響 . . . . . . . . . . . . . . . . . . . . . 32 3.3 3.4 第4章 4.3 4.4 シミュレーション 4.2.1 写り込み . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.2.2 にじみ出し . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 適応マスク法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.3.1 窓処理と適応マスク . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.3.2 適応マスクの原理と数式化 . . . . . . . . . . . . . . . . . . . . . 37 適応マスクを用いたブラインド・デコンボリューション . . . . . . . . . 38 4.4.1 画像情報のサンプリング . . . . . . . . . . . . . . . . . . . . . . 38 4.4.2 シミュレーション . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.4.3 サポート領域に関する考察 . . . . . . . . . . . . . . . . . . . . . 55 4.4.4 実験 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 まとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 シフトバリアント劣化像の復元 61 5.1 シフトバリアントな劣化 . . . . . . . . . . . . . . . . . . . . . . . . . . 61 5.2 シフトバリアント劣化像のブラインド回復 . . . . . . . . . . . . . . . . 63 4.5 第5章 5.2.1 連続分割法を用いた回復 . . . . . . . . . . . . . . . . . . . . . . 63 v 5.2.2 シミュレーション . . . . . . . . . . . . . . . . . . . . . . . . . . 65 5.2.3 実験結果 (一眼レフカメラ) . . . . . . . . . . . . . . . . . . . . . 72 5.2.4 実験結果 (携帯電話付属カメラ) . . . . . . . . . . . . . . . . . . . 75 まとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 動画の復元 79 6.1 フレーム処理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 6.2 実験 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.3 第6章 6.3 第7章 6.2.1 実験方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 6.2.2 PSF の推定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 まとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 結論 91 謝辞 95 参考文献 97 研究業績 103 1 第1章 序論 近年,携帯電話付属のカメラモジュールやカプセル型内視鏡に見られるように撮像素子 の小型化・高性能化の技術の向上が望まれている.これらの技術は撮像素子のハイテク化 に従い,撮影された画像の処理系等も複雑になってきており,高度な画像処理技術を必要 とすることが多い.しかしながら,光学系と画像処理系を同時に改善していくことは困難 な問題であり,今日盛んに研究が行われている分野ではあるが,画像処理系は光学系への 依存性が高く,デバイスごとに異なった画像処理系が必要となることが少なくない. 本論文では,画像回復においてロバスト性の高いブラインド・デコンボリューションの 手法を用いることにより,様々な光学系においてアルゴリズムを替える事無く画像を回復 するアルゴリズムを提案する.画像の劣化は一般的にはコンボリューション演算で表すこ とが出来る.ブラインド・デコンボリューション法とは,コンボリューション演算により 得られた劣化画像情報のみから物体の回復と劣化関数の推定を同時に行う逆問題として画 像処理をとらえた手法であり,元来の画像処理法では行う事が困難である劣化関数の情報 が得られない状況下において物体を回復する手法である.提案する手法の利点としては, 系を選ばずに回復が可能であると言うことから,無理に複雑な光学系を用いて画像回復に 適した画像を撮影したりすることなく,ごく簡単なモジュールで撮影された劣化画像の回 復を行うことが可能である.簡易なモジュールから撮影された劣化画像の回復を行えると 2 第 1 章 序論 言うことは,特に先に挙げた内視鏡や携帯電話のカメラモジュールを現在よりも小型化す ることが可能となり,さらなる精密化の発展を伺えることとなる. ブラインド・デコンボリューションを行うための基礎として,トランケートされた画像 に関しての考察が必要となってくる.本研究ではトランケートされた画像に対する前処理 として,適応マスクという手法を用いている.本手法は,画像が切り出された時に周りか ら写り込む光の影響を少なくするための手法である.一見するとあまり必要性の無さそう に見えるが,一般的に撮影されているスナップショットやビデオカメラの画像というもの はそもそもが周りからの光の写り込みが存在するものであり,その影響を無視することが 決して出来ないものである.現在までに行われているブラインド・デコンボリューション 法のシミュレーション法としては,周りが暗い物体を用いて劣化画像を作成している.ま た,劣化画像を考える際に実際に物体と劣化関数のコンボリューション演算から得られた 劣化画像を用いているために画像の縁部分に物体から漏れ出した光の情報をそのまま残し ている場合や,フーリエ面においてコンボリューション演算を行うことによりサポート領 域内で完璧なコンボリューション情報を持っている劣化像を用いている.また,実際の実 験においては,暗い実験室で撮影された画像や,天体画像などそれに準ずる劣化画像を用 いることによりトランケートによる情報の汚染や漏れを考慮に入れない系を考えているも のがほとんどである.一方で通常の劣化画像を用いて実験をしている研究に関してはその 結果としては決して芳しいものが得られているとは言えないものである.すなわち,ブラ インド・デコンボリューション法を用いた画像処理としては画像の縁,つまりトランケー トされた部分の情報を考えることが非常に重要となり,トランケートすることにより起こ る汚染・漏れを考慮に入れて処理を行う適応マスクは一般的な画像回復を行うにも必須の 技術とこの先なって来る可能性が高い. マスク処理としては窓関数をはじめとする形の変わらないマスクに関しては研究が行わ れ,信号処理の分野では頻繁に用いられている.しかしながら,画像処理に用いるための マスク処理はほとんどが矩形窓を用いたものであり,いわゆる窓関数というものも有用な ものが存在しない.適応マスク法は信号処理における窓関数を現在の状況により変化させ 3 る手法であり,推定された劣化関数から作成されるため様々な光学系に対して適用するこ とが可能である. 本論文は,適応マスクを用いたブラインド画像回復を提案し,その応用技術として,位 置により劣化関数が変化してしまうシフトバリアントな光学系により撮影された画像の回 復,および CCD カメラモジュールで撮影された動画の回復を試みる. また,トランケートされた劣化画像の回復が行えるという観点から,応用技術として物 体の位置により劣化関数が変化する光学系に適応マスクを用いることにより回復を行うこ とが可能となる.情報量が多く画像全体に対して回復を行うことが困難である動画に対し ても,切り取られた画像の一部から劣化関数を推定することが可能であるため,動画への 応用も例としてあげられる. 第 2 章では,基本的な画像処理の手法の説明をし,第 3 章ではブラインド・デコンボ リューションの基本的な最適化法の例を挙げる.第 4 章において,本研究の肝となる動的 適応マスクについて詳しく述べ,第 5 章・第 6 章においては,シフトバリアントな光学系 で撮影された劣化画像の回復と動的適応マスクの適用例として動画への応用について記述 する.第 7 章では,動的適応マスクを用いた画像処理のシミュレーション・実験結果をま とめ,本研究の総括を行う. 5 第2章 原理 2.1 画像処理 人間の得られる情報の中で 80%∼90% を占めているのが視覚情報と言われている.画 像処理は,視覚情報で得られる情報の 1 つである画像を様々な状況に応じて処置を施す 手法である.デジタルカメラを例に挙げると液晶モニターに映し出される画像はエッジ強 調処理や明るさ・ヒストグラムの補正を施した画像であるし,オートフォーカスにおける エッジ検出処理,顔画像認証のために行う前処理やノイズの低減処理,画像を保存する際 の圧縮処理と言った様々な箇所でその場に応じて全く違った過程から成り立つ画像処理が 行われている. これらのデジタル画像処理における画像は,圧縮方式によるファイル自体の記述に違い はあるものの,主に 2 種類の方法で表現される.1 つ目の手法は画像のピクセルをそのま ま行列の要素として表現する手法である.多くの場合,加法三原色 RGB を基準とした 3 要素に分けて表現される.画像の 1 ピクセルに対して,三原色各々の強度情報をそれぞれ 指定する手法である.2 つ目の手法は画像そのものではなく画像の持つ空間周波数で表し たものである.この手法は,座標 (x, y) の各色要素の強度分布に対して,(u, v) を空間周 6 第 2 章 原理 波数とおいたとき, F (u, v) = F[f (x, y)](u, v) Z ∞ = f (x, y)e−iux−ivy dxdy, (2.1) −∞ という変換を行ったものであり,この変換をフーリエ変換 (Fourier transform) と呼ぶ [1, 2].*1 物体空間からフーリエ変換を行い,扱う変数が空間周波数に変換されることによ り画像そのものの形ではなく画像の性質を知る事が可能となる.また,式 2.1 からわかる とおり,複素数となってしまうため関数の値としては物体空間における値に比べると扱い が複雑になる反面,コンボリューション定理 (重畳定理) と呼ばれる式変形が可能である ため,デジタル画像処理における重要性は非常に高いものとなっている.コンボリュー ション定理とは,簡単のため 1 変数で表すと, Z Z f (x0 )h(x − x0 )dx0 e−iωx dx F [f ∗ h] = Z Z 0 0 f (x0 )e−iωx × h(x − x0 )e−iω(x−x ) dx0 d(x − x0 ) Z Z 0 0 −iωx0 0 = f (x )e dx h(x − x0 )e−iω(x−x ) d(x − x0 ) = = F (ω)H(ω), (2.2) という関係式である.2 変数の場合においても全く同様に式変形することが可能である. 画像の表現手法は大体上記に挙げたとおりである.本研究における画像処理は画像処理 の中でも画像回復・画像再構成と言われるものであり,これは見えづらくなった劣化画像 を人間の目により見やすくするため,または機械的に画像認識などの処理を行いやすくす るために回復画像を求める手法である.現在アプリケーションとして使われることが多い アルゴリズムとして,ルーシー [3]・リチャードソン [4] 法がある. 一般的に劣化画像 g(x, y) は次の式の様に物体 f (x, y) と劣化関数 (以下 PSF:Point *1 1 次のフーリエ変換は通信における信号処理や音響解析などに用いられるものであり,この時は変数とし て ω = 2πν と行った周波数や角振動数が用いられる. 2.1 画像処理 7 spread function と記述する)h(x, y) の畳み込み積分の形で表すことが出来る. g(x, y) = h(x, y) ∗ f (x, y) + n(x, y), Z = h(x0 , y 0 )f (x − x0 , y − y 0 )dx0 dy 0 + n(x, y). (2.3) 演算子 ∗ はコンボリューションを,n(x, y) は加算ノイズを表す.式 2.3 の両辺をフーリ エ変換すると,式 2.2 より G(u, v) = H(u, v)F (u, v) + N (u, v), (2.4) となる.ここで,G, H, F, N はそれぞれ劣化画像,劣化関数,物体,加算ノイズをフーリ エ変換したもの,(u, v) は空間周波数を表す.一般的に画像処理と呼ばれるものは式 2.3 または式 2.4 で表現されている行列 g または G にフィルター処理を施すものである [5]. 先ほど挙げた物体空間における劣化像が式 2.3,周波数空間における劣化像が式 2.4 で表 されるものである. 例として,微分フィルターであるラプラシアンフィルターやソーベルフィルターと言っ たフィルター [6] の多くは式 2.3 の両辺に 3 × 3 もしくは 5 × 5 の簡単な行列を畳み込み 積分することによりエッジの画像を得るフィルターである.対して,ノイズ除去として広 く用いられている積分フィルターであるガウシアンフィルターやローパスフィルター等の 処理は式 2.4 で表されている周波数空間で行われることが多い.これら画像表現手法の違 いは,計算のコスト面,表現のわかりやすさ等の理由から使いわけがされているためであ り,本質的な違いは存在しない. 上記に挙げたフィルター例は,ごく簡単なフィルターである.ラプラシアンフィルター やソーベルフィルターは形には多少相違はあるものの,関数としてはどんな状況で用いる ものにしてもほぼ同じ関数の形になる.また,ガウシアンフィルターはサポートサイズが その大きな相違点になるが,パラメーターとしてはガウシアン分布の分散のみでそのほと んどの特性が決まってしまうフィルターであり,ノイズ減少のフィルターや,ボケフィル ターとしては非常に簡単に表すことが可能で,フーリエ変換しても関数の形が変わらない ことから多く用いられる.しかしながらこれらのフィルターはロバストな画像回復手法と 8 第 2 章 原理 しては用いる事が困難である. 2.2 デコンボリューションフィルター そこで考えられたのがインバースフィルターやウィーナーフィルターと呼ばれるデコン ボリューションフィルター [7] である.まず式 2.4 においてノイズが無視できる状況を考 える.このとき以下の様に式変形することが出来る. F (u, v) = *2 式 G(u, v) . H(u, v) (2.5) 2.5 で表されるフィルター処理がインバースフィルターと呼ばれる処理である.PSF が事前に測定できている状況や,PSF の推測が容易に可能であるときに用いられるフィ ルター処理であり,簡易な形ながらも良い精度での回復を行うことが可能である.その特 性から,画像回復を行いたい光学系と同じ光学系を用いて事前に各位置での PSF を測定 してしまうフィルターバンク手法に用いられることが多い.また,H(u, v) = 0 となって しまう成分の事を考慮し,少ないながらもノイズレベルを設定し, F (u, v) = G(u, v) , H(u, v) + α (2.6) と,ノイズパラメーター α(¿ 1) を分母に入れてしまう疑似インバースフィルターも広く 用いられているものである. 劣化画像の S/N が比較的低い場合には,こちらのインバースフィルターを用いること が不可能であり,インバースフィルターにノイズの影響の項を加えたフィルターがウィー ナーフィルターと呼ばれるものである.式 2.4 を以下のように変形する. · ¸ N (u, v) G(u, v) = H(u, v) + F (u, v), F (u, v) G(u, v) F (u, v) = . H(u, v) + Γ(u, v) *2 (2.7) (2.8) この表現方法であると,行列同士の除算になり正確に定義することができないため,本来は, F (u, v) = G(u, v)H ∗ (u, v) |H(u, v)|2 と記述するべきであるが,本論文中では理解のしやすさを考慮し,式 2.5 のような記述法を用いる. 2.3 ブラインド・デコンボリューション 9 ここで, Γ(u, v) = N (u, v) , F (u, v) (2.9) である.Γ というノイズの項を考慮に入れることにより,式 2.6 と式 2.8 の違いが大きく 出てきている分,ノイズの効果も同時に考えなくてはならなくなってくる.ノイズ効果の 項としてはベイズ理論を用いた統計的な分布を仮定し,劣化画像から求められたものを 使うのが広く行われている手法となっている.光学系の情報が得られる場合においては, ウィーナーフィルターが最も効果的かつ簡単な手法のひとつである. しかしながら,インバースフィルター・ウィーナーフィルターどちらにしても事前に光 学系の情報または劣化関数の情報が必要となってくる.一般的に撮像光学系を事前に知る ことは困難であり,たとえ光学系の情報が手元に合ったとしても,温度変化による空気の 揺らぎ等により光学系が少しでも変化してしまうと上記のようなデコンボリューション フィルターはそのままの形で用いることが不可能になってしまうという欠点もある. そこで考えられたのがブラインド・デコンボリューション法である. 2.3 ブラインド・デコンボリューション ブラインド・デコンボリューション法とは,式 2.3 に示されている中で,劣化画像 g(x, y) の情報のみから物体を回復させようと言う方法である.式の形を見ればわかるとおり,い わゆる不良設定問題であり,簡単に解くことは不可能である.そこで必要になってくるの が様々な最適化法であるが,こちらの紹介は次の章で行う.図 2.1 にブラインド・デコン ボリューション法の大まかなイメージを示した. 本節では,ブラインド・デコンボリューションが何故必要になってくるのかを述べる. まず,ブラインド・デコンボリューションの効果的な箇所としてあげられるのが,劣化画 像の情報とアプリオリな情報のみから物体の回復が行えるという点である.一般的に撮影 された劣化画像情報というものは,PSF の情報や,劣化した系の情報を直接知ることが困 難であるため,それらの情報が必要でないブラインド・デコンボリューション法が如何に 10 第 2 章 原理 Convolution PSF Object Blind deconvolution Degraded image 図 2.1 画像のコンボリューションとブラインド・デコンボリューションの関係. 効果的な手法であるかが伺える.アプリオリな情報として必要なものは単純なもので, 1. 画像を構成している成分は強度分布に依るものであるので,負値を持たない. 2. 物体・PSF は有限の空間内に存在する. という二点のみである.すなわち,必要になる情報は劣化画像のみとなり,他の情報が一 切不必要であることがわかる.言い換えると,極簡単なレンズと撮像素子だけの光学系は もちろん,ダブレットで構成されている単純な天体望遠鏡から複雑なレンズ系を組み込ん でいる一眼レフカメラ,実験室で撮影された CCD の画像から何気なく撮影されたスナッ プショットまでを同じアルゴリズムを用いて画像回復処理を行うことが可能であることが ブラインド・デコンボリューションの優れたところである. 一見,情報が少なすぎて不可能であると思われるが,劣化画像が物体と PSF のコンボ リューションで表すことが出来るのであれば,ゼロ・シートの考え方を用いることにより 解が唯一の組として存在することを Lane と Bates が数学的に検証・証明している [5, 8]. 2.4 まとめ Lane らは z 変換*3 した空間内で 2 次元の行列が唯一の形で因数分解されるという数学的 事実を用いてデコンボリューションを成功させた. Lane らによってコンボリューションによって合成された劣化画像が唯一の物体と PSF の組に復元できることが示されているので,ブラインド・デコンボリューションの有効性 が示されたと言うことになる. 2.4 まとめ 基本的な画像の表現手法とフィルター手法を挙げた. フィルター手法としては,物体空間におけるフィルタリング・周波数空間におけるフィ ルタリング,コンボリューションフィルターとデコンボリューションフィルターがあり, 各々場合に応じて用いられる事を述べた. またデコンボリューションは劣化関数が既知である場合に非常に有用な手法である事を 確かめた.デコンボリューション法の一種にブラインド・デコンボリューション法があ り,劣化関数 PSF が未知である時に用いられる,劣化画像情報とアプリオリな情報から 得られる拘束条件から回復物体を求める手法である. *3 ベクトルを多項式で表す方法の一つである.画像である 2 次元の行列は 2 変数の多項式に変換される. 例えば,n × n ピクセルのサイズを持つ画像は,2 元 n2 次方程式として表される. 11 13 第3章 最適化法 この章ではブラインド・デコンボリューションの最適化法についていくつか記述する. 3.1 最適化法とは 巡回セールスマン問題 (Traveling salesman problem: TSP) と呼ばれる NP 困難な問 題がある.この問題は,n 都市が存在し,どの 2 都市間の距離も定義されているとき,あ る都市からスタートして全ての都市を周り終えるまでに移動する距離を最小にする,とい う問題である.実際に全ての場合を列挙して解くことも不可能ではないが,都市の数が増 えるに従って実時間で解くことが困難になってくる問題である.さらに,「全ての場合を 列挙する」と言うことは明らかに解にはなり得ない場合の経路も考える事になってきて, コストは膨大なものとなる.TSP において厳密な答えを求めようと思えば,ある種のア ルゴリズムを用いて求めることは可能である.しかしながら,厳密な解ではなくとも経路 がなるべく少なくなるようにする解を求めたいときや,時間的コストをかけられないと き,もしくは実時間内で解く事が不可能であるときには,別のアルゴリズムが必要となっ てくる. そこで用いられるのが最適化法と呼ばれるアルゴリズムである.言い換えると,厳密な 解を解く必要が無い,もしくは厳密な解を求めるのが困難である時に用いられる手法であ 14 第 3 章 最適化法 り,いわゆるヒューリスティックアルゴリズムと呼ばれるものである.例えば,一元方程 式の解を求める際に使われるニュートン法なども最適化アルゴリズムの一つと言える.こ こで,いくつか最適化アルゴリズムの例を挙げていく. 3.1.1 最急降下法 Steepest descent method 評価関数の傾きを求め,変数を変化させることにより一番評価関数の値が小さくなる方 向,すなわち評価関数の勾配 (grad) が一番小さくなる変位を与えることにより,評価関 数を減少させていく方法.ある地点にボールをおいて転がし,ボールが止まったところを 最小値と見る手法.計算におけるコストが関数の勾配を求めることのみであるので,解を 得るまでの時間が非常に短い.一方で,単純に勾配の最小値だけを取ることになるので極 値の周辺からシミュレーションを始めてしまうと必ず極小値に陥ってしまうため,局所解 に陥る可能性が非常に高い. 3.1.2 山登り法 Hill climbing method : HC 山登りをするときに自分の周囲で常に一番高くなる方向へ進む事をモデル化した手法.現 在の解の近傍のから評価関数が一番小さくなるところを次の解とする方法.局所探索法の 最も単純なものとして知られる.こちらの手法に関しても,極値において必ず収束してし まうため局所解に陥る可能性が高い. 3.1.3 共役勾配法 Conjugate gradient method : CG[9, 10] 最急降下法を改良した手法.単純にポテンシャルが小さくなる方向ではなく,現在までの 修正ベクトルと共役になる方向へ変数を変化させる手法.最急降下法よりも計算の効率が 3.1 最適化法とは 高く,連立一次方程式の解法として良く使用されるもの.正確な解に収束するためにはあ る程度の前処理を必要とするものがあり,前処理付き共役勾配法と呼ばれる.解が存在す るときには確実にベクトルの次元回の計算を繰り返すことにより解に収束することが出 来る. 3.1.4 最尤法 Maximum likelihood estimation : MLE[11] 最尤推定法とも呼ばれる.ベイズ理論に乗っ取り,”最も尤もらしい値”を,与えられた データが従う確率分布の母数として求める手法.20 世紀初期の頃に統計学者ロナルド・ フィッシャーにより考案された.元々統計学の考えを基にしているため,上記に挙げたア ルゴリズムと違い正確な答えを得る事には決して向いているアルゴリズムではないが,よ り複雑な系を考える必要がある時には非常に有効なアルゴリズムである.非常に古いアル ゴリズムでありながら,今なお画像処理の分野では非常に良く用いられている優れた手法 である. 3.1.5 遺伝アルゴリズム Genetic algorithm : GA[12, 13] 自然界の動植物の遺伝子を解の候補と見立てて,遺伝子がより良い状況へと進化するのに 模倣しているアルゴリズムで,求めているものの大まかな関数を知るためには非常に高速 なアルゴリズム.解の候補として,固体を複数個用意し,適応度の高い固体を優先的に選 択,交差,突然変異の操作を繰り返すことにより解を探索する.選択・交差により遺伝子 の情報を改良していくが,選択・交差を続けると全ての遺伝子の情報が似たような情報と なってしまい,局所解に陥ってしまう.そこで,遺伝子に突然変異を起こすことにより, 局所解に陥る可能性を低くしている.比較的高速で局所解に陥りにくいアルゴリズムであ り,進化アルゴリズムの根幹を為すアルゴリズムとして広く研究が行われている. 15 16 第 3 章 最適化法 3.1.6 ニューラルネットワーク Neural network : NN[13, 14] 生物の脳機能に見られる特性をシミュレーションによって表現することを目指した数学モ デルで,シナプスの結合によりネットワークを形成した神経細胞が学習によってシナプス の結合強度を変化させる様子をモデル化した物である.ニューラルネットワークは大きく 分けて 2 つの手法があり,教師あり学習と教師なし学習がある.パターン認識などの目的 とするものが事前にわかっている場合は教師あり学習によりシナプスとニューロンが形成 される.一方でデータの分類やデータマイニングには教師なし学習が用いられる.計算量 は比較的少なく,良好な解を得られる事が多いため今日では様々なニューラルネットワー クの手法が確立されている.盛んに研究がされ始めたのは 1980 年代後半と最近のアルゴ リズムだと思われがちではあるが,比較的歴史は古く 1960 年代から研究されていた痕跡 も残っている. また,手法としては教師の有り・無しを別として大きく 3 つに分けられる. ・フィードフォワードニューラルネットワーク 一番初期に考案されたニューラルネットワーク手法で,入力層→中間層→出力層と言った 様に単方向への刺激が与えられるニューラルネットワークであり,画像認識や認証の分野 において広く用いられている.自己組織化マップ [15, 16] やパーセプトロン [17] がこの 手法に含まれる. ・フィードバックニューラルネットワーク フィードフォワードと違い双方向への刺激が存在するニューラルネットワーク.全ての ノードが結合しているものを全結合フィードバックニューラルネットワークと呼ぶ.ホッ プフィールド型ニューラルネットワーク [18] が代表的なものであり,画像回復手法に用い られることが多い. 3.2 ブラインド・デコンボリューションの最適化法 ・確率的ニューラルネット モンテカルロ法 [19] などのような確率的な刺激を与えたニューラルネットワーク.下記 に述べる焼き鈍し法もその 1 つと言える.ボルツマンマシン [20] がその代表例. 3.1.7 焼き鈍し法 Simulated annealing algorithm : SA[21] 金属を加工する際に,熱した後で徐々に冷やすことにより結晶を成長させ,中の不純物を 取り除いたり,結晶構造を整然なものにして欠陥を減らす作業のことを焼き鈍しという. この成長過程をモデリングした物が SA 法である.解を求めるための繰り返し計算の中で 温度パラメーター T を導入し,変数のグローバルな動きを可能とした,局所解に陥る状 況を作らないアルゴリズム.最初に T を高く設定しておくことにより局所解から抜け出 す事を容易にする.そして徐々に T を下げて行くことにより,解に収束していく事にな る.具体的な手法としては次の節で実際に挙げるが,上記 NN の中の考え方であるボルツ マンマシンや,メトロポリス・モンテカルロ法の考え方に基づいたアルゴリズムである. SA 法がその性質から非常に時間的なコストが高いにも関わらず画像処理や光学と言った 研究分野だけではなく,より広い研究分野において応用されている理由は,温度パラメー ターの減少速度を無限小,メトロポリス・モンテカルロ課程を繰り返す回数を無限大にす ることにより正確に厳密解にたどり着くことが証明されているからであり,これらのパラ メーターを緩和することによるエラーの増幅も他の手法と比べて小さいためである. 3.2 ブラインド・デコンボリューションの最適化法 さて,実際にアルゴリズムをいくつか挙げてきたが,中からブラインド・デコンボリュー ションに実際に適用したものの中から,第4章以降のシミュレーション・実験において実 際に用いているアルゴリズムについて紹介する. 17 18 第 3 章 最適化法 3.2.1 フーリエ反復法 (AD アルゴリズム) フーリエ反復法は位相回復法 [22] の一つとして良く知られる手法である.位相回復手 法における代表例は,エラーリダクション法や,ハイブリッドインプットアウトプット法 などがあり,これらの手法は物体空間と周波数空間でそれぞれ補正を行っていくという方 法である. g Fourier transform Satisfy Fourier constraints Satisfy function constraints g’ 図 3.1 |G|eiϕ Inverse Fourier transform |G’ |eiϕ 位相回復におけるエラーリダクション法.物空間と周波数空間におけるそれぞ れの拘束条件を別々にかけることにより,弱位相物体の位相を求める手法.強度情報を 物空間,位相情報を周波数空間においてそれぞれ拘束条件を与えることにより位相を回 復する. 例としてエラーリダクション法のブロック図を図 3.1 に示す.推定像 g をフーリエ変換 し,得られた G の強度情報へサポート外領域についての拘束条件を加え,G0 に更新する. その上で G0 を逆フーリエ変換することによりフーリエ空間で拘束条件を与えた物体空間 での像を得ることが出来る.得られた推定像は負値を含むため,負値をなくす or 軽減す る拘束条件を物体空間で与え,次の推定像 g を得る.そして次のイタレーションへと移る 作業を繰り返すことにより,G における正しい位相 φ を得るという手法である. 3.2 ブラインド・デコンボリューションの最適化法 19 フーリエ反復法をブラインド・デコンボリューションへ取り入れたのが Ayers と Dainty である [23, 24].Ayers と Dainty はフーリエ反復法においての物体空間での非負拘束条 件のかけやすさ,フーリエ空間におけるデコンボリューションの容易さに着目し,ブライ ンド・デコンボリューションへの応用を行った. Fourier transform F’ =G/H h Satisfy non-negative constraint for h f’ Feedforward h’ Inverse Fourier transform Inverse Fourier transform Satisfy non-negative constraint for f f H’ =G/F Fourier transform 図 3.2 Ayers と Dainty によって提唱されたフーリエ反復法を用いたブラインド・デ コンボリューション法 (AD 法).物体 f と PSFh の非負拘束条件を物空間で別々にか けることにより物体および PSF を求める手法. まず,入力された劣化像 g に対し,適当な PSF の初期推定像 h0 を入力する.h0 はフー リエ変換され,式 2.5 により物体のフーリエ像 F 0 を得る.F 0 は逆フーリエ変換され,物 体空間での推定物体 f 0 が求められる.ここで得られた f 0 は本来物体像が持つはずのない 負値を情報として持っているため,負値を減らす or 無くすような拘束条件をかけること により,負値が減少した or 存在しない推定像 f を得る.h に関しても同じような処理を 施し,フーリエ変換,デコンボリューション,逆フーリエ変換,拘束条件を適用と繰り返 すことにより,推定 PSF h を得る.この一連の操作を繰り返すのが AD アルゴリズムと 呼ばれている手法である.途中に存在するフィードフォワードパラメーターは,f および h の補正が前のイタレーションにおける関数とあまりに大きくかわらなくするようにする 20 第 3 章 最適化法 ための補正値であり,F 0 と F ,H 0 と H をそれぞれある割合で混合するためのパラメー ターである. 説明を見てわかるとおり,AD アルゴリズムはフーリエ変換と逆フーリエ変換,フーリ エ空間内での単純な割り算,そして極簡単な拘束条件を適用すると言う非常に単純な計算 の繰り返しであるので,その計算コストは低く,非常に高速で解を得られるアルゴリズム になっている.その反面,解の収束性は保証されているものではないので,解の信頼性に 関しては高いとは言えない.しかしながら,その高速性は他の多々あるブラインド・デコ ンボリューションのアルゴリズムの中では群を抜くものになっており,広く使用されてい るアルゴリズムでもある. 本手法は 2 人の頭文字をとって一般的に AD アルゴリズムと呼ばれている. また,フーリエ反復法を用いたブラインド・デコンボリューション法は他にも研究され ており,先ほど例にだしたエラーリダクション法とハイブリッドインプットアウトプット 法を用い,画像の強度分布とフーリエ空間内での位相を同時に扱う様な違う観点から用い たブラインド・デコンボリューション法も存在する [25]. 3.2.2 ゼロシート法 ゼロシート法 [5, 8, 26] とは,z 変換を基にしたブラインド・デコンボリューション法 である.他の手法が逆問題として最適化法を用いて近似解を求めるのに対し,ゼロシート 法は唯一確実な解を求める手法として知られている. z 変換 [27] とは,信号処理に用いられる変換法の 1 つであり,f (x, y) を画像情報とし たときに F (u, v) = ∞ X ∞ X f (xm , yn )u−m v −n , (3.1) m=0 n=0 とする変換である.ここで (xm , yn ) は (x, y) を離散化表現したものである.形を見てわ かるとおり,x, y のサンプリング間隔を t としたとき x = xm = mt, y = yn = nt となる 3.2 ブラインド・デコンボリューションの最適化法 21 ため,u = eiξt , v = eiηt と置換する事により t → 0 としたときにフーリエ変換そのもの の式となる. さて,画像の性質として,サポート領域が有限 (N × N ) である事と画像情報が 2 次元 であるが必要となってくる.このとき,画像の z 変換は F (u, v) = N X N X f (x, y)ux v y , (3.2) x=1 y=1 と書くことが出来る.このとき劣化像 g の z 変換 G は, G(u, v) = F (u, v) · H(u, v), (3.3) と表される.ここで 2 変数の多項式が実数体上において因数分解出来るときには唯一の 2 つの多項式の積に分離することが出来ると言う数学的性質を用いてブラインド・デコンボ リューションを行うのがゼロシート法である.実際の過程としては G(u, v) = 0 になる点 を 4 次元空間内でプロットしていき,プロットした曲面がなめらかな 2 つの閉曲面とな るように分離をしていく.4 次元空間を用いると言うことで計算が非常に複雑になり,ま た物体のサポート領域が広がるに従って指数関数的に考えなくてはいけない変数が増加し てしまうために計算コストが非常に高くなってしまう弱点がある.さらに,ノイズが乗っ てしまうと本来なめらかである閉曲面同士の交わっている部分においてなめらかでは無く なってしまうために分離が難しくなってしまうと言う短所もある.ただし,それをさしお いても正確な解を得られるという点から,ゼロシート以上のブラインド・デコンボリュー ション法は存在しないと言っても良い. また,位相回復の分野へゼロシート法を応用する研究もなされている [28]. 3.2.3 焼き鈍し法 (SA 法) 前節で述べた SA 法をブラインド・デコンボリューションに適用する [29].まず SA 法 の説明で述べられていたメトロポリス・モンテカルロ法について説明をする. 22 第 3 章 最適化法 モンテカルロ法 [19] とは,物質中の中性子の運動を測定するためにジョン・フォン・ノ イマンにより考案された手法であり,数値計算する段階において乱数を用いて行う手法で ある.例えば,ある領域 (Region of interest : ROI) の面積を求めたいときにシンプソン 法や台形近似,三角近似などを良く行う.これに対してモンテカルロ法では,ROI を丸ご と含んでしまうような面積 S である領域を準備し,その中に無数の点を乱数により n 点 打点する.打点された中から,ROI に含まれた点の数 m を測定する.n が十分に大きい ときには,乱数が一様であるならば確率的に ROI の面積は mS/n になるというのがモン テカルロ法による面積の求め方である.感覚的にもわかるとおり,n が大きくなるほど測 定精度は向上する. 一方メトロポリス法 [30] とは,モンテカルロ法によって得られた結果を受け入れる際の 確率を考慮する手法である.系のエネルギーを E とし,モンテカルロ法により E が dE だけ変化するとしたときに, dE ≤ 0, (3.4) である (エネルギーが減少した) 状態であればモンテカルロシミュレーションによる数値 の変化を 100% の確率で受け入る.もし, dE > 0, (3.5) µ ¶ dE p = exp − , kB T (3.6) である (エネルギーが増加した) 時は, と言う確率に従いモンテカルロシミュレーションの数値の変化を受け入れるような手法を 言う.ここで,T は系の温度,kB はボルツマン定数を表す.また,ここで現れた受け入 れ確率にボルツマン分布を用いている事から,一連の流れをボルツマンマシンと呼ぶ.例 えば,ホップフィールド型ニューラルネットワークにボルツマンマシンを組み込むことに より,学習がスムーズに行くことが知られている. ボルツマンマシンにおいて dE < 0 における受け入れ確率 p は温度で左右されることに なるが,最初に系のエネルギーと温度を高く設定しておき,徐々に温度を下げることによ 3.2 ブラインド・デコンボリューションの最適化法 23 り金属における結晶をきれいに並べ不純物を取り除く過程に見立てた数値解析法をシミュ レーテッドアニーリング法 (焼き鈍し法) と呼ぶ. ここで,SA 法をブラインド・デコンボリューションに適用する [29] ことを考えてみる. 以下にブロック図を示す. 䊝䊮䊁䉦䊦䊨ᴺ䈮䉋䉍 fࠍᄌൻ 䊝䊮䊁䉦䊦䊨ᴺ䈮䉋䉍 hࠍᄌൻ g’ =f㧖h g’ =f㧖h E=||g’ 㨺g|| E=||g’ 㨺g|| 䊜䊃䊨䊘䊥䉴ᴺ䈮䉋䉍 fࠍᦝᣂ 䊜䊃䊨䊘䊥䉴ᴺ䈮䉋䉍 hࠍᦝᣂ ᷷ᐲᦝᣂ no no 䊗䊦䉿䊙䊮䊙䉲䊮䈱࿁ᢙ䈏 ⷙቯ࿁ᢙ䈮㆐䈚䈢䋿 yes ᷷ᐲᦝᣂ࿁ᢙ䈏 ⷙቯ࿁ᢙ䈮㆐䈚䈢䋿 㫆㫉 䉣䊈䊦䉩䊷୯䈏 චಽ䈮ਅ䈏䈦䈢䋿 yes ജ䋧ᓳర⚳ੌ 図 3.3 焼き鈍し法を用いたブラインド・デコンボリューション法.E は最適化の評価 関数である.図の中では単純なノルムとなっているが,系のエネルギーと増減が同じに なり得るものであればどんな形でもかまわない. まず,初期予想物体と初期予想 PSF を用いて g 0 をコンボリューションにより作成し, エネルギーの初期値を E = kg 0 − gk として,画像復元を開始する.初期予想物体 f に対 してモンテカルロシミュレーションを行うことにより f の関数を少し変化させる.具体 24 第 3 章 最適化法 的には,f の乱数によってピクセルを抽出する.抽出されたピクセルの値を乱数により少 し変化させる.このとき,変化量も温度関数依存にすることにより,温度パラメーターが 高い時点に収束速度を上げることが可能である.次に,変化した f と,前のイタレーショ ンから求められている h をコンボリューションすることにより,現段階での予想劣化像 g 0 ,さらにそれを基にエネルギー E を求める.現エネルギー値 E と前段階のエネルギー 値から,エネルギーの差 dE を求め,それを基にメトロポリス法に従い f の更新をする. PSF h に対しても同様の作業を繰り返す. h の更新が終了したら,ボルツマンマシンの回数が同一温度の状態で規定回数繰り返さ れたかどうかを判定し,繰り返し回数が足りないようであれば温度を変えずにもう一度 f の更新から繰り返す.ボルツマンマシン過程を繰り返した回数が十分であるのならば,次 に温度の更新回数およびエネルギー値の確認を行う.温度の更新回数が規定回数に達して いなかったり,エネルギー値が終了判定エネルギー値よりも高かった場合には,温度更新 を行った上でボルツマンマシン過程に戻る.温度の更新回数が規定回数に達するか,エネ ルギー値が十分小さいと見なされたときには,解析を終了し結果を出力する. 以上が SA 法を用いたブラインド・デコンボリューションの手順となる. 上記では f ,h を更新する際にランダムにピクセルを抽出して行うと記述したが,解析 の初期のうちはピクセルの抽出はせずに全てのピクセルに対してボルツマンマシン過程を 繰り返した方が収束の速度は上がる. AD アルゴリズムと比べ,1ピクセルごとに作業を行う必要があるため非常に解析に時 間がかかるが1ピクセルごとに少量ずつ変化を繰り返して最適化を行っているので安定し た収束性と高精度な推定を得ることができる. 3.2.4 ハイブリッド AD/SA 法 そこで,上記に挙げた二つのアルゴリズム,AD 法と SA 法をハイブリッドしたアルゴ リズムを考える.AD アルゴリズムは解の収束性が保証されていない一方で,非常に高速 3.3 シミュレーション に近似解を得ることが出来ることから,物体の大体の形を得るために計算の初期に導入し ている.一方で,SA 法を用いるのは,大体の物体の形を得ることが出来た後,より正確 な解を得るために用いている.これは前節においても説明したが,SA 法は解の収束性が パラメーターによりある程度保証されているためである. この手法は,SA 法の高速にある程度の復元を得られる利点と SA 法のより精度の高い 復元を得られる利点を併せ持ったアルゴリズムである.手順としては,AD アルゴリズム を最初に劣化画像に適用し,まずは大体の物体と PSF の関数を推定する.その後,SA 法 の初期予想として AD アルゴリズムで得られた物体と PSF を代入することにより,高速 かつ安定した収束性を得ることの出来るアルゴリズムを考えることが出来る. 本論文では 4 章以降ではブラインド・デコンボリューション法として,ハイブリッド AD/SA 法を用いることとする. 3.2.5 その他のアルゴリズム 他のアルゴリズムの例としては,現在 MATLAB の Image processing toolbox でも採 用されている方法である最尤法 [31, 32],ウェーブレット変換を用いたもの [33, 34],ホッ プフィールドニューラルネットワーク [35, 36, 37],遺伝アルゴリズム [38, 39],ハイブ リッド GA/SA[40],遺伝アルゴリズムを改良した免疫アルゴリズムよる最適化 [41] など, 多岐にわたる研究が活発になされている.最近ではイタレーションのない最適化法 [42] も考案されているが,こちらは解析時間はほぼ単純なフィルター処理のみなので非常に高 速ではあるが,得られた回復像には擬像や負値が多く含まれていることもあり,これから のさらなる研究を望まれる分野である. 3.3 シミュレーション ここで実際にブラインド・デコンボリューション法を用いて画像を回復する.前章に挙 げられたうち,ホップフィールドニューラルネットワークを基にしたボルツマンマシンを 25 26 第 3 章 最適化法 用いたブラインド・デコンボリューション法を例としてシミュレーション結果を示す. 3.3.1 原理 ホップフィールドニューラルネットワーク (HNN) は相互結合型のニューラルネット ワークであり,ユニットの更新によりネットワークが平衡状態に達するアルゴリズムで ある. i ニューロンから j ニューロンへの結合定数を Tij とすると,Tij = Tji 6= 0 を満たす ネットワークである.i ニューロンの出力は他のニューロンからの出力 νj の総和とバイ アス入力 Ii から,k + 1 回目のイタレーションの出力 νi (k + 1) は, input(i) = n X Tji νj (k) + Ii , (3.7) i (input(i) > 0) 1 νi (k + 1) = νi (k) (input(i) = 0) 0 (input(i) < 0) , (3.8) となる.ここで n はニューロンの数とする.HNN の各状態におけるエネルギーは, n n n X 1 XX Tij νi νj − Ii νi , E=− 2 i j i (3.9) により表現される. 劣化像 g ,物体 f をベクトル表示し,PSF を畳み込み行列 H により表現したとき,劣 化像を基準とした現在の推定像・PSF のエネルギーは以下の式で表すことが出来る. E= 1 2 kHf − gk . 2 (3.10) 物体のサポート領域を N × N として式 3.10 を計算すると, 2 E= N 1X 2 gl − l=1 2 2 N X hl,m fm , (3.11) m=1 となる.ただし,hl,i は物体 fi が劣化像 gl に与えるインパルス応答とする.ここで fi = n X k=1 νi,k , (3.12) 3.3 シミュレーション 27 Convolution 図 3.4 ブラインド・デコンボリューションのシミュレーションに用いる画像.左側の 2 つをコンボリューションすることにより右側の劣化像を作成し,求められた劣化像か ら左側の 2 つの画像を得ることを目的とする. と表されることを考えると,式 3.9 と式 3.12 を比較することにより, Tij = − Ii = X X hl,i hl,j , (3.13) l gl hl,i , (3.14) l とすれば良いことがわかる. HNN では基本的に発火するかしないかだけで考えられてしまうので,出力として νi (k + 1) = 1 , 1 + exp(−input(i)) (3.15) のシグモイド関数を用いてこの値を画素の値として出力する.また,HNN ではエネル ギー値が減少する方へしか更新は受け入れられないため,局所解に陥る可能性が出てきて しまう.そこで,更新時にはメトロポリス法を用いてある程度のエネルギー上昇を受け入 れるようにする.物体の推定が終わったら PSF の推定へ移る.物体と PSF の推定を何 回か繰り返すことによりシミュレーションを行う. 28 第 3 章 最適化法 図 3.5 ボルツマンマシンを用いて行ったブラインド・デコンボリューション法のシ ミュレーション結果.図 3.4 の劣化像情報から,2 つのコンボリューションの要素を求 めた. 図 3.6 ボルツマンマシンを用いて行ったブラインド・デコンボリューション法のシ ミュレーション結果.図 3.4 の劣化像に 15% の加算ガウシアンノイズを乗せた劣化像 から,2 つのコンボリューションの要素を求めた. 3.3.2 結果 シミュレーションに用いた物体,PSF を図 3.4 に示す.ボルツマンマシンを用いて行っ たシミュレーションの結果が図 3.5 である.また,15% のガウシアンノイズを乗せた劣化 像から得られた結果を図 3.6 に示す.ノイズが乗っているため,回復像にもノイズが乗っ てしまっているが,元の図形がくっきりと見えるのがわかる. 3.4 まとめ 3.4 まとめ 一般的な不良設定問題の最適化法から画像処理における最適化法のうち代表的なものを 紹介した.画像回復アルゴリズムに用いられる最適化法からいくつか例をあげ,ブライン ド・デコンボリューション法への応用として一般的に用いられているアルゴリズムについ て,その手法と考察を第 2 節においてまとめた.また第 3 節では,実際のブラインド・デ コンボリューション法として第 2 節で例として挙げたボルツマンマシンの手法を用いてシ ミュレーションを行った. 第 3 節の結果から非常に複雑な不良設定問題であるブラインド・デコンボリューション 法を用いた劣化画像回復の可能性と,その有効性を実証した. 29 31 第4章 適応マスク手法 4.1 従来手法のブラインド・デコンボリューションの問題点 前章では様々な最適化法とブラインド・デコンボリューションにおける代表的な最適化 法を挙げ,多くの参考文献を挙げた.しかし,これらの文献で用いられている劣化画像に 関して一つ重大な欠点が残されているのである.実際に撮影された画像として多くの論文 で挙げているのは天体画像である.元々の研究で天体画像のスペックル画像の画像復元を 目指していたのでこれに関しては致し方ないことであると思われる.次に多いのが暗室の 中で撮影された画像である.前章のシミュレーションを行った画像については周りが黒で 埋められている画像である. ここに共通している欠点が何であるかというと,周りからの光の写り込みの影響をほと んど考えなくても良い状況下であると言うことである.シミュレーションを行っている 研究に関しても同様で,前もって撮影された画像に対して PSF を仮定してコンボリュー ションを行い作成された劣化像を用いて解析を試みている.こちらに関しても画像の縁の 部分は何の処理を施すこともなく,基の画像からにじみ出ている光をそのままの状況で残 し,サポートの外側は真っ暗な状況を意識することなく想定している.しかしながら,一 般的にこのような撮影像が存在するかというと,天体画像や特殊な画像以外にはほとんど 無いと言っても過言ではない. 32 第 4 章 適応マスク手法 そこで我々がブラインド・デコンボリューションを行う際に考慮しなければ行けない課 題は,周りからの写り込みのある劣化画像を回復する事の出来る処理と言うことになる. 例えば,画像を分割したときに各々の画像を回復する事を考えてみる. 従来,そのような研究が全くされていないわけではなく,例えば画像を分割してそれぞ れを並行してニューラルネットワークにかけることにより画像の復元をしている研究 [43] 等がある.ここに挙げた参考論文を実際に見てみると,とてもではないがうまく画像が復 元出来ているとは言い難く,分割した画像の縁部分を全く考慮に入れていないため,復元 画像に非常に多くの擬像が現れてしまっている.また,画像をトランケートすることに より限られた周波数領域におけるブラインド・デコンボリューションを行っている研究 [44, 45] もある.こちらの研究では,セグメント化を行いお互いに隣同士のセグメントに 与える影響を考え非常に良い結果が得られているが,セグメント化したにもかかわらずそ の一番の利点である計算コストはあまり変わらない結果となっている. このように,トランケートされた劣化画像に関する画像復元法は必要な技術にも関わら ずあまり研究がなされておらず,解決法がほとんど無い状態になっている.また,逆に言 うとこの研究分野が解決することにより,画像処理技術の様々な面において応用をするこ とが出来,画像処理技術の発展には今後欠かせない問題になってくるとも考えられる. 4.2 トランケートが画像に及ぼす影響 そこで,本論文では,周りからの写り込みの影響を考えたブラインド・デコンボリュー ション法を提案する.実際には,物体と PSF のコンボリューションにより作成された劣 化画像をトランケートすることにより,周りのピクセルから PSF によって写り込みの影 響が大きい劣化画像セグメントを抽出し,抽出された劣化画像を回復することにより計算 機シミュレーションを行う.また,CCD カメラにより撮影されたテストチャート画像か らトランケートすることによりセグメントを抽出し,抽出されたセグメントから PSF を 求めることにより物体全体を復元する実験を行う. 4.2 トランケートが画像に及ぼす影響 33 4.2.1 写り込み 実際にトランケートが画像に及ぼす影響を挙げていく.まず第一点目が周りの物体から の光の写り込みである. 光の写り込みの様子を図 4.1 に図示した.まず,物体に関して考える.左図の黒く塗り つぶされた範囲がトランケートされた画像領域 (ROI),その周りの白い領域がセグメント に隣接する ROI 外の領域である.次に,劣化画像について考えると,ROI と ROI 外の 領域が別々に劣化する訳ではなく,同時に劣化していることは式 2.3 によりわかる.右図 でもわかるとおり,白い部分と黒い部分が劣化により混ざってしまっている.このとき, 元々黒く塗りつぶされた領域に重なってきているグレーの領域が外からの光の写り込みに よる影響の部分である. 一般的に撮影された画像に対しては,必ず縁の部分にこのような汚染が存在するはずで あり,デコンボリューションを行う際には決して無視する事が出来ない効果を及ぼしてい Object Degraded image Degrading Trancated region (ROI) Neighboring pixels of ROI Contamination from neighboring pixels 図 4.1 PSF による劣化が原因で起こる ROI 周辺のピクセルから ROI への光の写り 込みによる影響.図中の白い領域から PSF による劣化により,黒い領域 (ROI) へ光 が写り込んでいる (グレーの領域). 34 第 4 章 適応マスク手法 る.従来提案されているフィルター処理やデコンボリューションでは,この”写り込み領 域”を考慮に入れられていないものがほとんどであり,それ故に特に画像の縁部分での劣 化が良く見られるのである. もちろん,劣化画像をトランケートやセグメンテーション化する際には,図からわかる ようにそこには多分の「余計な情報」が含まれていて,この余計な情報を除去できない限 りセグメント内の劣化画像のみからのデコンボリューションによる画像復元は困難である ことがわかる. 4.2.2 にじみ出し 先ほどの写り込みでは「周囲からの余計な情報」が問題であった.一方,逆に考えると 「周囲への情報の拡散」も避けられない問題であり,考慮に入れる必要がある. 図の見方は先ほどと同じである.ROI からの情報の拡散が見て取れるであろう. こちらの情報に関しては前記した周辺ピクセルからの余分な情報の写り込みの影響ほど Object Degraded image Degrading Trancated region (ROI) Neighboring pixels of ROI Leak to neighboring pixels 図 4.2 ROI から周辺ピクセルへの情報のにじみ出し.図 4.1 とは逆に,今度は ROI(白 い領域) から周辺のピクセル (黒い領域) への情報の漏れ (グレーになっている部分) が 多く存在している. 4.3 適応マスク法 大きな影響を及ぼすわけではないが,やはり復元画像の縁部分に関してはうまく処理して やらないと行けない部分である.こちらに関しては,物体のサポート領域を狭めてやるこ とによりその影響を軽減させることが出来る. 先にも述べたとおり今日発表されている研究成果のほとんどが天体画像や実験室内で周 りをなるべく暗い状態にして実験されているものがほとんどである.このような画像に対 しては従来の手法を何の前処理をすることなく適用しても周りからの光の写り込みが極端 に少ないためにシミュレーションと同様に画像復元をすることが可能となっている.ま た,上記二つの情報による汚染・拡散を考慮に入れないとトランケートされた画像のデコ ンボリューションはうまくいかないことは,参考文献 [43] に載せられている復元画像を見 ればわかるうえ,ここで重要になってくる点は,ごく普通に撮影された画像は全てが一般 的に撮影の過程においてトランケートされているという点である.本論文においても,こ の章の最後の節において,実際に情報の汚染・拡散を考慮に入れていないときの復元画像 と適応マスクにより特に情報の汚染を軽減した時の復元画像を比較する. 4.3 適応マスク法 前節で挙げられた二つの大きな劣化画像の復元における問題要因を軽減させる方法を考 える. 4.3.1 窓処理と適応マスク トランケートされた信号の処理系統として,窓処理 [46] というものが存在する.窓処理 とは,トランケートされた信号に関してダイナミックレンジと周波数特性を欲しいデータ に合わせて変化させるような一種のマスキング処理の事である.通常,窓関数は一次元の デジタル信号処理に用いられることが多く,信号の用途によりハミング窓やハニング窓と 言った様々な窓関数が適用される. 35 36 第 4 章 適応マスク手法 窓処理とは具体的には,窓関数と言うものを用いて定義される.窓関数からいくつか代 表的なものを挙げてみる.まず,矩形窓と呼ばれるものは w を窓関数とし, w(x) = 1, (0 ≤ x ≤ 1), (4.1) で表される.着目する領域内を 1,それ以外の領域の情報を無くしてしまうと言う非常に 簡単なものである.これに対してよく使われるハニング窓,ハミング窓はそれぞれ次の関 数を用いて表現される. w(x) = 0.5 − 0.5 cos 2πx, (4.2) w(x) = 0.54 − 0.46 cos 2πx. (4.3) 式 4.2 がハニング窓,式 4.3 がハミング窓と呼ばれる窓関数である.また,まとめて w(x) = a − (1 − a) cos 2πx, (4.4) とする一般化ハミング窓と称される窓関数も存在する.ただし,0.5 ≤ a ≤ 1 で定義され るものとする.式 4.4 は a = 1 で矩形窓,a = 0.54 でハミング窓,a = 0.5 でハニング窓 となる. 実際に行われる窓処理とは,信号 f (x) に対して上記に挙げられた窓関数 w(x) を, h(x) = w(x) × f (x), (4.5) と掛け合わせ,得られた h(x) の周波数特性を見ることにより信号の特性を測定する手法 の事を言う. 画像処理に関しては現在のところ特殊なケースをのぞいて利用されることは少ないが, 基本的に画像と呼ばれるものは意識すること無く矩形窓を用いて表現されている.矩形窓 の特性上,周波数解像度は理論的に全ての窓関数において最も高くなっているが,ダイナ ミックレンジに関しては低いものとなっている. 本論文では,窓処理を画像処理に適用し,さらに一般的に用いられているハミング窓や ハニング窓など一定のパラメーターのみで表される窓関数ではなく,状況に応じて動的に 4.3 適応マスク法 37 変化するようなマスキング処理を行うことにより,前節で挙げられた周辺物体からの光の 写り込みを減少する事を目的とする.この窓処理に用いられるマスクを適応マスク [47] と呼ぶ. 4.3.2 適応マスクの原理と数式化 それでは実際に適応マスクがどのような働きをするかの説明とマスクの数式化について 述べる. 本研究の中心となる技術である適応マスクは,周辺ピクセルからの光の写り込みの影響 を少なくすることを目的として作成される.周辺物体からの光の写り込みを減らすために は,画像の縁部分に周りからの PSF の写り込みを減らすようなマスク処理をしてやれば よい.そこで本研究においては,推定された PSF を用いてマスクを作成し,それを劣化 画像のセグメントにかける事を提案する.マスクを作成する時点での PSF を h とおい て,マスク m(x, y) は以下の式に従い作成される. m(x, y) = d(x, y) ∗ h(x, y). (4.6) ここで,d(x, y) は物体と同じサポート内で定義されたダミーオブジェクトであり,一様 関数となっている.つまり,ダミーオブジェクトと PSF の畳み込みの形によりマスクを 作成する.この形の物理的な意味としては,すなわち,物体のサポート領域内から PSF によってにじみ出る領域をグレースケールのレベルまで考慮した形で考える,と言うこと であり,サポート領域内から PSF によって劣化した画像がセグメントのサポート領域内 のどこまで影響を与えているかを考えている.逆に言うと外側からの写り込みをなるべく 多くシャットアウトするような形のマスクを作成しており,周辺ピクセルからの光の写り 込みによる汚染を減少する役割を果たしている. そして作成されたマスクは劣化画像のセグメント情報 s(x, y) に対して, sm (x, y) = s(x, y) × m(x, y), (4.7) 38 第 4 章 適応マスク手法 という形で適用される.sm (x, y) は,マスク処理されたセグメント情報であり,sm (x, y) を実際にブラインド・デコンボリューションのアルゴリズムに代入することとなる. 4.4 適応マスクを用いたブラインド・デコンボリューション 適応マスクを用いたブラインド・デコンボリューションの効果をシミュレーションと実 際に CCD カメラで撮影された画像に関して実験を行う. 4.4.1 画像情報のサンプリング 実際にシミュレーション・実験を行う前に,ここまでで考えられている関数は全て連続 関数で考えられているため,画像のサンプリングをしなくてはならないのでそれについて 述べたいと思う. まず最初に,式 2.3 は離散化したときに, g[k, l] = (f ∗ h)[k, l] = FN X FN X f [m, n]h[k − m, l − n], (4.8) m=1 n=1 h[m, n] = 0, if [m, n] ∈ / F N 2, (4.9) と表される.[k, l] は離散化されたピクセルの座標であり,物体のサポート領域を F N × F N の矩形領域として計算を行っている.同様に本研究で用いる疑似インバースフィル ターは, F [s, t] = *1 と表すことが出来る.α *1 1 G[s, t], H[s, t] + α (4.10) はノイズレベルであり,ノイズレベル未満の値を無視する役割 離散化した場合も第2章*2 での脚注と同様,厳密に言うと本来は F [s, t] = G [s, t] H ∗ [s, t] |H [s, t]|2 + β と式に表すべきであるが,簡単のために式 4.10 の様な表現方法をとる.離散化した場合では,各行列の 要素ごとに割り算をすると考えてば良い. 4.4 適応マスクを用いたブラインド・デコンボリューション 39 をしていると同時に,フィルターの分母が 0 になることを防ぐパラメーターでもある. 一方サンプリングされたマスク処理はどうなるかというと, mask[k, l] = (d ∗ h)[k, l] = FN X FN X d[m, n]h[k − m, l − n], m=1 n=1 1 d[m, n] = FN2 1 1 ··· 1 1 1 ··· .. .. . . . . . 1 .. . 1 1 ··· 1 (4.11) , (4.12) sm [k, l] = mask[k, l] × s[k, l] (4.13) と表すことが出来る.ただし,実際のシミュレーションを行う際には,mask 関数ががこ のままだとダイナミックレンジが非常に低くなってしまっているので, mask[k, l] = (d ∗ h)[k, l] + (1 − max {(d ∗ h)[k, l]}) = FN X FN X d[m, n]h[k − m, l − n] m=1 n=1 Ã + 1 − max ( FN FN XX )! d[m, n]h[k − m, l − n] , (4.14) m=1 n=1 と補正することにより, ¡ ¤ mask[k, l] ∈ 1 − max(d ∗ h), 1 (4.15) として,マスク処理された劣化画像セグメント sm [k, l] のダイナミックレンジの減少を防 ぐ.図 4.3 に実際に作成されたマスクの例を示す. 4.4.2 シミュレーション 以上で述べた適応マスクの効果を確かめるためにシミュレーションを行う. この小節では,物体と PSF をそれぞれ用意し,それらをコンボリューションすること により劣化画像を作成し,作成された劣化画像から一部分をトランケートすることにより 40 第 4 章 適応マスク手法 Currently estimated PSF Dummy object 図 4.3 Adaptive mask 適応マスクの例.適応マスクは式 4.12 で表されるダミーオブジェクトとマス クを求める時点での推定 PSF から作成される.初期マスクは PSF を平均値フィル ターとして計算をする. セグメント情報を得る.得られたセグメント情報に対して初期適応マスクを適用し,初 期マスク処理されたセグメント情報に”適応マスクを組み込んだ SA 法に基づくブライン ド・デコンボリューション法”を用いて数値計算を行う.エネルギー関数 E は,推定 PSF h・復元物体 f をコンボリューションした物およびマスク処理されたセグメント情報 sm より以下の式により求める. E = k(h ∗ f )[k, l] − sm [k, l]k sX n £ ¤2 o (h ∗ f )[k, l] − sm [k, l] = (4.16) k,l 本研究においては,物体そのものの情報ではなく PSF を求めることを目的とする.何 故ならば PSF を求めることが出来れば,もとの劣化画像に対して式 4.10 の疑似インバー スフィルター処理を行うことにより全体の復元像を得る事が可能であるからである. ここで重要になってくるのがマスクの更新処理になる.マスクの更新は SA 法の温度 処理に従い行う.SA 法の温度パラメーターを 10 回更新するたびにその時点で推定され 4.4 適応マスクを用いたブラインド・デコンボリューション ている PSF h を用いてマスクの更新を行う.マスクの更新を行ったら,次のメトロポリ ス・モンテカルロ処理に入る前に劣化画像のセグメント情報 s に対してマスク処理を行 い,新しいマスクでマスク処理されたセグメント情報 sm をメトロポリス・モンテカルロ 過程に代入する.毎 10 回の温度更新の時点でマスクは新しいものへと更新される.以上 のシミュレーション過程をブロック図で表したのが図 4.4 である.SA 法を用いたブライ ンド・デコンボリューションに,マスク更新処理とその判定過程を入れたものである. 以上でシミュレーションの準備が整ったので,実際にシミュレーションを行う.用いた 物体は,Lenna *2 ( 図 4.5) を用いる.256 × 256 ピクセルのサイズで 8bit のグレースケー ࡑࠬࠢᦝᣂ ࠣࡔࡦ࠻ᖱႎᦝᣂ yes ࡔ࠻ࡠࡐࠬ ࡕࡦ࠹ࠞ࡞ࡠㆊ⒟ ᷷ᐲᦝᣂ࿁ᢙ䈏 ⷙቯ࿁ᢙ䈮㆐䈚䈢䋿 㫆㫉 䉣䊈䊦䉩䊷୯䈏 චಽ䈮ਅ䈏䈦䈢䋿 no ࡑࠬࠢᦝᣂ್ቯ ᷷ᐲᦝᣂ no yes ജ䋧ᓳర⚳ੌ 図 4.4 適応マスクを用いたシミュレーテッドアニーリング法に基づくブラインド・デ コンボリューション法のブロック図.マスクは温度更新が 10 回行われるごとに更新さ れる.図のマスク更新判定は温度更新が 10 の倍数になるかどうかと言う分岐である. *2 画像処理シミュレーションにおいて頻繁に使われている画像である.肌の部分のなめらかな領域,帽子の 41 42 第 4 章 適応マスク手法 図 4.5 Lenna.適応マスクを用いたブラインド・デコンボリューション法のシミュ レーションに用いる物体.256 × 256 ピクセルの大きさである. ルで表されている. PSF として,適当な標準偏差を持たせたサイズ 7 × 7 ピクセルの正規分布の画像 (図 4.6) を用いる.ただし,上記正規分布は小数点以下第3桁目を切り捨てる (丸める) 形に して,わざと丸め誤差を持たせている.完全なガウス関数の形にしないのは,ガウス関数 同士の積がガウス関数になってしまう,つまりガウス関数が何回コンボリューションをし てもガウス関数になってしまう関数であり,ブラインド・デコンボリューションの解とし て厳密な形でのガウス関数が入ってしまうと解が無数に存在してしまうからである.小数 羽に見られる複雑な高周波成分をもった領域,そして髪の毛の様に一方向のみに高周波を持つ領域と言っ た画像処理を行う際に最も必要とされる要素を兼ね備えた画像である故,多くの画像処理の研究に用いら れている. 4.4 適応マスクを用いたブラインド・デコンボリューション 点以下第3桁目と言うとほとんど誤差の出ない感じを受けるかもしれないが,グレース ケールを 8bit で表しているので,第3桁目の切り捨てと言うと無視出来ない誤差を持っ ている.そのため,ガウス関数とはかなりの違いが出てくる. 図 4.6 シミュレーションに用いた PSF.7 × 7 ピクセル. 図 4.5 および図 4.6 に図示された物体と PSF をコンボリューションすることにより劣 化画像 (図 4.7) を作成し,この劣化画像を回復することにより適応マスクの効果をシミュ レートする.まずはじめに,図 4.7 から,20 × 20 の領域をトランケートする.図 4.8 に 示されているように,図の中心から少し右にずれたあたり,Lenna の”左目から鼻の頭に かかる部分”をセグメントとして抽出した. 初期マスクは,PSF としてアベレージフィルター [6, 7] を用いダミー物体とコンボ リューションすることにより求めた図 4.9 を用いる.マスク処理されたセグメント sm は 図 4.10 に示すように,セグメントの縁の部分がかなり暗くなる様な図になっている. 43 44 第 4 章 適応マスク手法 図 4.7 本シミュレーションで用いる劣化像.図 4.5 と図 4.6 の画像をコンボリュー ションしたもの.この劣化像を回復する事を目的としている. 図 4.8 図 4.7 から切り取られた 20 × 20 ピクセルのセグメント.このセグメントのみ を用いて PSF を推定する. 4.4 適応マスクを用いたブラインド・デコンボリューション 図 4.9 シミュレーションに用いた初期マスク.初期マスクはアベレージフィルターと ダミーオブジェクトのコンボリューションにより作成される.本シミュレーションに おいては,7 × 7 のアベレージフィルターと 14 × 14 のダミーオブジェクトから作成さ れた. 図 4.10 初期マスク処理を施されたセグメント.セグメント (図 4.8) と初期マスク (図 4.9) の要素同士の積. 45 46 第 4 章 適応マスク手法 図 4.11 適応マスクを用いたシミュレーションにより求められた物体. 図 4.12 適応マスクを用いたシミュレーションにより推定された PSF. 図 4.10 に示されたマスク処理された劣化像を SA 法に代入することにより,復元され たトランケートされた物体が図 4.11 に,同様に推定された PSF が図 4.12 に示されてい る.ここでは,2000 回の温度更新を行う事を終了条件としてシミュレーションを行った. 図 4.6 と図 4.12 を比べてみることにより本手法の効果が非常に高いことがうかがえるで あろう. 参考として,最終的なマスキングされた劣化セグメントを図 4.13 に示す. 先ほども述べたとおり,本研究で求めたいものは PSF の情報である.そこで PSF がど れだけ回復できているかを知る定量的な指標として,二乗平均誤差 (Mean square error: MSE) を導入する.MSE は式 4.17 と表される. MSEh = HN HN 1 XX 2 (h[m, n] − h0 [m, n]) HN 2 m=1 n=1 (4.17) ここで,PSF のサポート領域を HN × HN とし,h0 [m, n] は図 4.6 そのものを表してい 4.4 適応マスクを用いたブラインド・デコンボリューション 図 4.13 47 最終的な適応マスク処理を施された劣化セグメント.図 4.11 と図 4.12 をコ ンボリューションした物とほぼ一致している. る.図 4.12 の PSF の MSE は, MSEh = 1.59 × 10−4 , (4.18) となっており,8bit のグレースケールで言うとピクセル全体を通して平均 7∼8 階調ほど のずれにとどまっている. さて,セグメント情報から PSF がかなりの精度で推定できる事がわかったところで, 本シミュレーションで行っている劣化が位置によって劣化関数が変化しない,いわゆるシ フトインバリアントな劣化過程であることに着目する.シフトインバリアントと言うこと は劣化画像全面にわたって PSF が同じ関数であると言うことであり,劣化画像のサポー ト領域内であればどこからセグメントを抽出しても同じ PSF が推定できるはずである. そこで,抽出するセグメントを増やしたシミュレーションを行ってみる. 先ほどと同じ領域を抽出したものと,先ほどとは別の,”口の右下あたり”をトランケー トしてセグメントとして抽出する (図 4.14).そしてこれらの2つのセグメントを同時に 解析することにより,PSF の推定の高速化・高精度化を試みる. まず,同じ処理をした上でシミュレーションを開始する.2つのセグメントが同時に解 析されていると言うことは2つの PSF の推定が得られていることになる.もちろん,10 48 第 4 章 適応マスク手法 図 4.14 もう 1 カ所のセグメント.Lenna の口のあたりをトランケートしたもの. 回の温度更新が終わった時点でも PSF の推定が2つ得られることになる.便宜上それら を PSF1 ( h1 ) と PSF2 ( h2 ) とおく.PSF1 と PSF2 から, ha [k, l] = 1 (h1 [k, l] + h2 [k, l]), 2 (4.19) として PSFa を計算する.これら3つの推定 PSF に対してエリート選択処理を行い2つ の PSF を残す.残された2つの PSF からマスクをそれぞれ作成し,ボルツマンマシン過 程に代入する. と言う処理を行った.この結果得られた回復物体が,図 4.15 である.同様に得られた 推定 PSF が図 4.16 である.このとき PSF の MSE は, MSEh = 7.35 × 10−5 , (4.20) となっており,先ほどの半分以下,グレースケールに直すとわずか平均 3 階調ほどの誤差 まで精度が上げられている. と,記述されてもよくわからないと思うので,次のページに比較として回復物体および 推定 PSF をそれぞれ並べてみる.図 4.17 は物体の復元像である.(a) がシングルセグメ ントの時,(b) がダブルセグメントでの回復像,(c) はマスク処理をせずに,さらに PSF を正解のものに固定したときに SA 法により求めた回復像,(d) は物体そのものを表して 4.4 適応マスクを用いたブラインド・デコンボリューション 図 4.15 セグメントを 2 つ用いたときの回復セグメント画像. 図 4.16 セグメントを 2 つ用いたときの推定 PSF. いる.図 4.17(b) と (d) を比べてみてわかるように,セグメントを2つ用いて回復された 像は非常に高精度であることがわかる.また (c) からわかるとおり,適応マスクを用いて いないときは,ほとんど回復されない.この結果から適応マスクの有効性が示された.ま た同様に推定された PSF も比較として図 4.18 に載せる.こちらは目で見て判別できるほ どの違いが無い事がわかると思われる. 49 50 第 4 章 適応マスク手法 図 4.17 (a) (b) (c) (d) セグメント画像の比較.(a) 1 つのセグメントから求められた回復物体.(b) セグメントを 2 つ使って求められた回復物体.(c) 適応マスクを用いずに劣化セグメン トのみから SA 法によって求められた回復物体.正しい PSF をシミュレーションに用 い,PSF の推定は行っていない.(d) 現物体のセグメント. 4.4 適応マスクを用いたブラインド・デコンボリューション (a) 図 4.18 (b) 51 (c) 推定 PSF の比較.(a) 1 つのセグメントから求められた推定 PSF.(b) セグ メントを 2 つ用いたときに求められた推定 PSF.(c) 現 PSF. 図 4.19 評価関数 E のグラフ.セグメントを 2 つ用いることにより,収束性が大幅に 向上している. 図 4.19 はエネルギーの値をセグメント1つを用いたときとセグメント2つを用いたと きそれぞれに関して記したものである.縦軸にエネルギー値,横軸には温度更新回数で表 している.セグメント1つを用いて復元したとき 2000 回の温度更新が必要であったレベ ルまで,セグメント2つを用いた場合は約 500 回で復元することが可能であった.また, 大体同じ計算回数になるシングルセグメントで温度 2000 回更新とダブルセグメントで温 52 第 4 章 適応マスク手法 度 1000 回更新を行った場合についてのエネルギー値を比べてみると,シングルの場合は 1.49,対してダブルの場合は 0.52 となっている.シングルとダブルで約 2.9 倍の差が出 てきている.これは同じ時間で考えるとダブルセグメントで考えた方がより優れた回復を 見せると言うことを実証している.このグラフからもダブルセグメントを用いる有用性が 実証された. 参考までに,図 4.20 は実際に人間の目に映る状況を考えて,各手法で得られた復元像 と実際の物体からトランケートしたセグメントを線形補間して並べたものである. (a) 図 4.20 (b) (c) 人間の目から実際見るときの状況を近似して回復物体を線形補間した画像. (a) はシングルセグメント,(b) はダブルセグメントにより求められた回復物体である. (c) は現物体セグメントを線形補間したもの. 最後にセグメントを2つ用いた時に求められた推定 PSF から疑似インバースフィル ターを作成し,作成されたフィルターから劣化像全体を回復したものを図 4.21 に掲載し ておく.全体像としても非常に良く復元できていることがわかる. 4.4 適応マスクを用いたブラインド・デコンボリューション 図 4.21 適応マスクを用いて得られた推定 PSF により逆フィルター処理を施して劣化 画像全体を回復させたもの. 53 54 第 4 章 適応マスク手法 図 4.22 PSF のサポートサイズを変化させたときのエネルギー関数のグラフ.HN = 7 が正しいサポートサイズであり,このときエネルギー関数が最小になっていることが わかる.注目する点は,HN ≤ 7 である時には収束値は減少し続けるが HN = 8 と なったときにエネルギー値が急激に増加する点である. 4.4 適応マスクを用いたブラインド・デコンボリューション 55 4.4.3 サポート領域に関する考察 実際にシミュレーションを行ったわけであるが,1つだけ大きな問題が存在する.先ほ どのシミュレーションを行う際には,物体と PSF のサポート領域を完全に固定して考え てしまっている.もちろんシミュレーションであるから可能な作業である. そこで,この問題点を解決するために,物体・PSF のサポートサイズを変化させエネル ギーを測定したものが図 4.22 である. サポートのサイズが間違ったサイズであると,エネルギー値が大きな値で収束している ことがわかる.特に,PSF のサポートサイズが少しでも大きい状態になるとエネルギー 値はかなり大きな値になってしまう.そのため,同じ状況下においてイタレーションを何 度か繰り返し,サポートサイズを変化させていく.変化させていく中で急激なエネルギー の上昇がある前後で正解のサポートサイズが存在することになる.また,サポートサイズ が現実のサポートサイズよりも小さい場合においてもエネルギーの上昇が見られた.これ らは 1 だけ替えた場合に限らず,さらに変化させていくと徐々にエネルギーが上昇して いったため,現実のサポートサイズで一番エネルギーが小さくなると考えて良い.ただ し,PSF のサポートサイズが 1 ピクセルになることだけは考えてはいけない.この状態 では,劣化画像がどのような画像であるかに関わらず,デルタ関数と劣化画像のコンボ リューション,すなわち, g[k, l] = δ ∗ g[k, l], (4.21) と,分解できてしまうためである.因数分解や素因数分解において,1 を例外にするのと 同じ理由である. この理由からサポートの領域を多少大きくしても自然とエネルギーが下がると PSF の サポート領域は決まってくるので,物体・PSF どちらのサポート領域も大きめに取ること により,サポートの評価をすることが無く解析を行えることになる. 56 第 4 章 適応マスク手法 4.4.4 実験 実際に撮影された画像の回復することを考える.用いるのは図 4.23 に示されたテスト チャートのピンぼけ画像である. この図からセグメントを抜き出し,実験を行った.実際に撮影された画像にはノイズの 効果が入ってしまうので,前処理として,5 × 5 ピクセルのメディアンフィルターを適用 した.コンボリューション線形フィルターを用いない理由はもちろん,劣化画像に別の関 数をコンボリューションしてしまうと,解をユニークに求めることが不可能になってしま うからである. 前小節の手法により推定した PSF を用いてインバースフィルターにより復元した像が, 図 4.24 に表されている図である.先ほどのシミュレーションの時に比べて図の情報量が 大きいので,SA 法だけで戻すことは時間的コストが非常に高くなってしまう.それ故, 本実験においては,最適化法に AD アルゴリズム [23] を取り入れている.AD アルゴリ 図 4.23 実験に用いる劣化像.テストチャートのピント外れ劣化像を用いる.撮影に 用いたのは,デジタル一眼レフカメラ (EOS-Kiss Digital). 4.4 適応マスクを用いたブラインド・デコンボリューション 57 ズム過程において用いられるエネルギー関数は SA 過程とは違うものであり,物体空間に おける負値成分の和の割合を考えたものとした. sP EAD = Pk,l k,l 2 [k, l] f− f 2 [k, l] , (4.22) ただし,f− は推定像 f の負の成分とする.この AD 過程の操作を入れることにより,大 まかな物体・PSF 両方の形を知ることにより,PSF の推定を実時間内に終わらすことが 可能となる. さて,どこまで細かいところまで復元できているかを見るために一部拡大図を載せる. 図 4.25 に示されているのはテストチャートの Line 3 の No. 6 の部分の画像である.マス クを使わずに求めたときの復元像とマスクを用いて求めたときの復元像を並べたものであ り,マスクを用いたときはラインがかなり高い精度で復元できていることが見て取れる. 図 4.24 図 4.23 の劣化像をハイブリッド AD/SA 法に基づいた適応マスクを用いて回 復した画像. 58 第 4 章 適応マスク手法 (a) (b) 図 4.25 劣化像をトランケートしたセグメントから求められた回復物体から抜き出し たセグメント.セグメントはテストチャートの 3 ブロック目の 6 番目の部分.(a) 適応 マスクを用いて回復したもの.(b) 適応マスクを用いずに回復したもの. 4.5 まとめ 適応マスク手法について述べ,AD/SA 法に基づいた適応マスクを用いたブラインド・ デコンボリューション法のシミュレーションおよび実験を行った.シミュレーションで 特に注目するべき箇所は,図 4.17(a),(c) の比較である.正しい PSF に固定してシミュ レーションを行ったにもかかわらずマスク処理を施していない場合はほとんど回復物体を 得ることが出来てない一方,適応マスクを用いることにより図 4.17(a) の様に物体が回復 していることがわかる.これらの図の比較により適応マスクの必要性とその有効性が示さ れた.また実験においてもテストチャートの画像を見てわかるとおり,非常に高精度な回 復が行われていることがわかる.図 4.25 の比較からも,適応マスクを用いたときは同じ 計算回数の中で高精度な画像を得ることが可能である. 適応マスクがブラインド・デコンボリューション法において必要不可欠な技術であるこ 4.5 まとめ とがシミュレーション・実験の双方から実証された. 59 61 第5章 シフトバリアント劣化像の復元 前章ではトランケートされた劣化像の回復を行った.本章と次章では,トランケートさ れた劣化像のブラインド回復の応用技術について触れたいと思う. 本章では,トランケートされた劣化像が回復出来るという利点から,画像を小さい領域 に切り分け,各劣化像に対してブラインド回復を行うことにより,位置によってインパル ス応答が変化してしまうような画像の回復を試みる. 5.1 シフトバリアントな劣化 シフトバリアントな劣化関数とは,物体の位置 (ピクセル) により劣化関数が変化して しまうような劣化関数のことを言う.例として,シフトバリアントな劣化関数による劣化 の様子を図示する.まず,シミュレーションで使う画像を図 5.1 とする.こちらの画像は 8 × 8 ピクセルのサイズの画像である. このとき,例えば中心の 4 × 4 ピクセルとそれ以外の外側のピクセルに対してそれぞれ 違う PSF をかけると言う考え方がシフトバリアントな劣化関数による劣化である.例え ば,それぞれの領域に対する劣化関数を次のように定義してみる. この場合,中心部分はそれほど劣化はしないが,外側の部分での劣化が大きい様な PSF に設定されている.なお,図 5.2 の中では違いをわかりやすくするために中心部分を 1 に 62 第 5 章 シフトバリアント劣化像の復元 (b) (a) 図 5.1 シフトバリアントな PSF による劣化の回復シミュレーションに用いる物体. (a) 図 5.2 (b) シフトバリアントな PSF.物体ピクセルの中から,中心領域 4 × 4 ピクセル部 分が (a) により,それ以外の外側の領域は (b) によって劣化する. 5.2 シフトバリアント劣化像のブラインド回復 63 規格化しているが,実際は光の強度エネルギーが 1 になるように規格化して計算を行う. 実際にシフトバリアントな劣化関数,図 5.2 により図 5.1 を劣化させて見たのが図 5.3 で ある.シフトインバリアントな劣化関数による劣化 (コンボリューションでの劣化) に比 べて劣化の仕方が複雑であることがわかる. (b) (a) 図 5.3 シフトバリアントな PSF 図 5.2 により劣化した劣化像.(a),(b) はそれぞれ 図 5.1(a),(b) の物体に対応している. 5.2 シフトバリアント劣化像のブラインド回復 5.2.1 連続分割法を用いた回復 前節で説明したシフトバリアント劣化像のブラインド回復を行う.数値解析を行う前 に,シフトバリアントな劣化関数による劣化をモデリングしておく.離散化した表現で書 くと, g[k, l] = XX m n f [k − m, l − n]hk,l [m, n], (5.1) 64 第 5 章 シフトバリアント劣化像の復元 となる.ここで,hk,l は f [k, l] に対するインパルス応答である.また,疑似コンボリュー ション行列を用いることにより, g = Hf = [h0 h1 · · · hn ]f , (5.2) と表すことも出来る.f ,g は物体,劣化像をそれぞれベクトル表示にしたもの,hi は物 体の i 番目の要素にかかるインパルス応答であり, hi→0 hi→1 .. . hi = hi→j .. . hi→N , (5.3) と表し,hi→j は物体 f の i 番目の要素から劣化像 g の j 番目の要素へのインパルス応答 と定義する. 現在までのところ,小出,小松によって特異値分解を用いて H の疑似逆行列を求める ことによりブラインド復元を行う行列反復演算法 (MIM) と言う研究が発表されている [48].MIM は n 枚の劣化像から, G = [g0 g1 · · · gn ] , (5.4) F = [f0 f1 · · · fn ] , (5.5) と置き,行列演算の中で非負拘束条件とサポート条件を用いた最適化を行い物体行列 F と PSF 行列 H を交互に推定していく手法である.疑似逆行列を用いることにより非常に 精度の高い回復をすることに成功している.しかしながら,この手法においてはその計算 の処理上,同じシフトバリアント系で劣化された画像が,劣化関数行列 H のランクと同 程度必要になる.例えば,先ほど挙げた図 5.3 のような劣化画像を回復するためには,36 個程度の劣化画像が必要となってくる.必要な劣化画像の数は,物体のサポートサイズが 上がれば上がるほど増大していき,計算機的なコストはそれほど大きなものとはならない 5.2 シフトバリアント劣化像のブラインド回復 が,同じ状態で撮影された劣化画像の数が同程度増えていくので,実用的な手法とは言え ない.もちろん,顕微鏡画像など何枚も同じ状況で撮影される事が保証されている光学系 においては非常に有用なものとなる. そこで,前章で紹介した適応マスクを用いた手法を用いることによりトランケートされ た劣化像を回復することにより,各々の領域でのインパルス応答を求めることにより物体 を復元することを目的として,連続分割法 [49, 50] という手法を提案する. 連続分割法は劣化画像の領域をフェーズに従って縦横半分ずつに,すなわち 1/4 に刻ん でいく手法である.刻んだ縁部分はもちろん周辺のピクセルからの光の写り込みが存在す るために,そのままの状態ではうまく回復することが困難であるので,ここで適応マスク を適用することとなる.詳しく書いたのが図 5.4 である. ただし,劣化画像の切断を行う場面において,各部分の PSF のサポート領域よりセグ メントが小さくなってしまうと回復することが出来なくなってくるので,PSF のサポー ト領域よりセグメントが十分大きいと見なせなくなったときには,式 5.2 の各成分に対し て SA 法を用いることにより物体の復元を行う手法に変更する. この連続分割法は,物体のあるピクセルに注目したとき,そのピクセルの周辺のピクセ ルに関しては劣化関数は似たものになる事に基づき,ある領域内ではインパルス応答はピ クセルによらず一定の関数であると見なせるという過程に基づいて行っている.そのた め,劣化関数が急激に変化する部分においては,ある程度のフェーズが進み,デコンボ リューションを行う領域がある程度小さくならないとうまく回復をすることが出来ず,注 意が必要となってくる. 5.2.2 シミュレーション 図 5.3 で示されている劣化画像の回復を実際に行ってみる.まず,全体的な流れをブ ロック図にして図 5.5 に表す. 最初に行われているゼロパディングとは,領域を 2n にするための作業である.この作 業は,2 分割ずつにしていくので領域を分けやすくすると共に,フーリエ変換をかけると 65 66 第 5 章 シフトバリアント劣化像の復元 2 f 1*h 1 f 2*h 2 f*h 2 f 3*h 3 0th phase 4 f 4*h 4 1st phase N f 12 f 11 *h 11 *h 12 4 f 14 f 13 *h 13 *h 14 2nd phase 図 5.4 ... N Nth phase 連続分割過程におけるフェーズ.各フェーズにおいて,領域を 1/4 ずつに抽出 する.各セグメント内での PSF はシフトインバリアントとし,異なるセグメント間で は PSF は別物として仮定する.連続分割は物体のサポート領域が 1 ピクセルになるま でか,エネルギー値が十分小さくなるまで繰り返される. 5.2 シフトバリアント劣化像のブラインド回復 67 Blurred image Zero padding Segmenting Iterative BD procedure Recovered object 図 5.5 連続分割法における最適化のブロック図. きに手間を省くために行われる.次に,適応マスクを用いたブラインド・デコンボリュー ションを行う.こちらに関しては前章で用いたものと同じようにイタレーションの初期で は AD 法を,その後に SA 法を適用する方法をとっている.そして各セグメントに一定回 の温度更新を行ったら,次のフェーズのために分割を行う.この作業を繰り返し行うこと により,画像の復元を行った.1フェーズ分を実際に図 5.3(a) を用いて表したものが図 5.6 になっている. 以上の操作を行った最終的な復元像が図 5.7 である.二乗平均誤差 (Root mean square error: RMSE) s Erms = |f − f0 |2 |f0 |2 (5.6) をモニタリングした物を図 5.8 に示した.また,最終的に得られたシフトバリアント PSF を図 5.9 に示す. 68 第 5 章 シフトバリアント劣化像の復元 2n 図 5.6 (a) (b) (c) (d) 連続分割法のフローチャートを図で表したもの.(a) 第 1 過程:サポートの大 きさが 2 の累乗になるまで周囲に値が 0 であるピクセルをパディングする.この図の 場合であると一辺が 24 = 16 ピクセルになるまでゼロパディングを行う.(b) 第 2 過 程:AD アルゴリズムをゼロパディングされた劣化像全体に対して適用する.(c) 第 3 過程:第 2 過程または第 4 過程において得られた物体およびゼロパディングされた劣 化像をそれぞれセグメント化する.(d) 各セグメントは各々 AD 法または SA 法を用い て最適化される.最適化により求められた物体は第 3 過程へと戻り,反復計算される. 5.2 シフトバリアント劣化像のブラインド回復 69 (a) 図 5.7 連続分割法によって回復されたシフトバリアントな PSF により劣化した劣化 像の回復像.(a),(b) はそれぞれ図 5.3(a),(b) に対応する. (b) 70 第 5 章 シフトバリアント劣化像の復元 RMS error 1.0 phase 2 phase 3 0.1 phase 2 0.04 0 図 5.8 10 20 30 40 time (s) 50 60 70 式 5.6 で表される RMSE をモニタリングしたもの.実線はハイブリッド AD/SA 法を用いたもので,AD 法は反復回数を 20 回を 1 フェーズ,SA 法は 200 回の 温度更新 2 フェーズを行っている.点線は SA 法のみで行ったシミュレーション結果. 5.2 シフトバリアント劣化像のブラインド回復 71 図 5.9 (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) (m) (n) (o) (p) 連続分割法により得られたシフトバリアント PSF.それぞれ物体の 8 × 8 ピク セルに対して,(a) 左上 2 × 2 ピクセルに対する PSF,と言うように 2 × 2 ピクセルを 単位として対応している.物体の内側の 4 × 4 ピクセルに対する PSF が (f),(g),(j), (k) であり図 5.2(a) に,他のものが外側の PSF 図 5.2(b) に対応している推定 PSF で ある. 72 第 5 章 シフトバリアント劣化像の復元 5.2.3 実験結果 (一眼レフカメラ) 前小節において,連続分割法を用いた有用性が実証された.そこで,実際にカメラで撮 影された画像が回復できるかどうかを確かめてみる.撮影に用いたカメラは Canon EOS Kiss digital lens kit である.ここでは,シフトバリアントな系の代表例として,撮像素 子から距離が大きく異なる物体の写真を撮影し,その画像の復元をすることを試みる.撮 影における光学系は図 5.10 に描かれているとおり,カメラから前面にある物体までの距 離が約 1m,後方においてある物体までの距離が約 3m になるように設定されている. 実際に撮影された画像と復元された画像を図 5.11,図 5.12 に示す.まず,撮影された 画像 5.11 を見ると,前方・後方に存在する物体の両方にピントが合っていないことがわ かる.この画像に連続分割法を用いて数値計算してやることにより図 5.12 のように,前 方・後方ともにピントの合った画像になっていることがわかる.ただし,適応マスクの実 験を行ったときと同様,ノイズが存在するために前処理として画像全体に 5 × 5 ピクセル のメディアンフィルター処理を施している.結果としては,前後の境目や,後方の壁のな 3m Second target Camera First 1m target 図 5.10 シフトバリアントな PSF による劣化画像を撮影するための光学系. 5.2 シフトバリアント劣化像のブラインド回復 図 5.11 図 5.10 に示された光学系においてデジタル一眼レフカメラによって撮影され た劣化像.画像全体の大きさは 1024 × 1024 ピクセルにトランケートされている. 図 5.12 シフトバリアントな PSF による劣化像 (5.11) をハイブリッド AD/SA 法に 基づいた適応マスクを用いた連続分割法によりブラインド回復した物体. めらかな部分などに多少擬像が生じているが,本手法が非常にうまく画像を復元出来るこ とが実証された.また,シフトバリアントな系として PSF が求められている事が目で見 てわかるように図 5.13 に前後それぞれの部分から一部分のインパルス応答を示す.ここ でインパルス応答は見やすさのために適当に規格化されていることに注意してほしい. 73 74 第 5 章 シフトバリアント劣化像の復元 図 5.13 連続分割法により得られた推定 PSF.代表的な位置として,手前・奥の物体 にかかるインパルス応答を 1 カ所ずつ挙げた. 5.2 シフトバリアント劣化像のブラインド回復 図 5.14 携帯電話 (NTT Docomo 社の SH505i) 付属のカメラを用いて撮影されたシ フトバリアント PSF による劣化像.撮影条件は一眼レフデジタルカメラで撮影したも のとほぼ同じような状況である. 5.2.4 実験結果 (携帯電話付属カメラ) 似たような光学系において携帯電話付属のカメラを用いた画像についても画像の回復を 試みた.撮影に用いられたのは NTTDocomo SH505i である.まず,撮影された画像を 図 5.14 に示す. そして復元された画像が図 5.15 である. 先ほどの復元物体に対して,携帯電話で撮影された画像に関してはかなりのノイズが 乗っていることがわかる.この理由としては,携帯電話の中から生データとして画像情報 を取り出すことが不可能であり,携帯電話の機能として,撮影された画像が JPEG イメー ジへ変換されてしまっているためである. JPEG イメージに変換されることにより,高周波成分は消える事になるのだが,逆にコ ントラストの高いところにおいて高周波成分が多く出てしまっている.この理由としては JPEG イメージのコーディングの仕方にあり,JPEG イメージはある程度の大きさ,例え ば 8 × 8 ピクセルごとに画像ファイルを切り分けその中で高周波成分を消すことにより画 像を圧縮している.それ故にコントラストが高い部分においては,JPEG イメージの中に 75 76 第 5 章 シフトバリアント劣化像の復元 図 5.15 携帯電話付属の CCD カメラモジュールで撮影された劣化像・図 5.14 を適応 マスクを用いた連続分割法により回復した物体. 図 5.16 図 5.14 の回復物体に対する PSF の分布.グレースケールが似た値になるほ ど近い PSF に,グレースケールが大幅に異なるところは PSF が大幅に異なっている 事を示している. 5.3 まとめ 77 は,擬像が現れたりしてしまう.また,ブロックノイズと呼ばれる JPEG や MPEG に存 在するノイズも現れる.こちらは,先ほどのピクセルに区切った時点で隣の領域との不整 合が出るためであり,画像全体としてみると望んでいない周波数成分を持つことが多く, この周波数成分によるノイズをブロックノイズと呼ぶ.さらに,携帯電話内での処理とし て,画像におけるエッジの先鋭化処理が施されているのが図 5.14 から読み取れる.これ らの原因から,図 5.15 における回復像は多くのノイズを含んでしまったり,エッジ先鋭 化フィルターによる擬像をデコンボリューションしているのでフィルターが原因で起こる エッジ部分の擬像を含んでしまったものだと考えられる. さて,携帯電話で撮影された画像だが,前方と後方の劣化関数がどれだけの範囲で散ら ばっているかを見るために図 5.16 を掲載する. v uP P ¡ ¢2 u k l hc (k, l) − h(m,n) (k, l) t Q(m, n) = 1 − P P 2 k l (hc (k, l)) (5.7) 式 5.7 は各位置 (m, n) における PSF h(m,n) と中心の位置における PSF hc の違いを表 したものであり,図 5.16 は Q の分布を表しており,PSF 自身が近いほど近いグレース ケールになるように描かれている.くっきりと前方と後方において劣化関数が変化してい ることがわかるであろう. 5.3 まとめ 適応マスク手法の応用として,物体の位置により劣化関数が変化してしまうシフトバリ アントな PSF による光学系により劣化された劣化像の回復法の提案をした. シフトバリアントな劣化系による劣化の回復を行うに当たり,比較的位置が近いところ 同士は PSF がそれほど変化しないと言う仮定の下,連続分割法という手法を考案した. 連続分割法は,劣化像・推定物体のそれぞれのサポート領域を 1/4 の領域に連続的に分割 していく手法であり,各セグメント内ではシフトインバリアントな系であると仮定してブ ラインド・デコンボリューションを行う手法である. まず,シミュレーション結果から適応マスクを用いた連続分割法によるシフトバリアン 78 第 5 章 シフトバリアント劣化像の復元 トな系の劣化像の回復が行えることを実証した. 実験には,一眼レフデジタルカメラで撮影された画像と携帯電話付属の CCD カメラを 用いて撮影された画像を用いた.どちらの場合も高精度な回復を行うことが出来たが,携 帯電話付属の CCD カメラに関しては,ファイルの保存をする際に JPEG 画像に変換さ れてしまっているために,多くのブロックノイズや不自然な縦の分割など見直すべき点が 多く残されている. 79 第6章 動画の復元 前章では適応マスクの応用としてシフトバリアントな系の劣化像の像回復を行った.そ れに対して本章ではアプローチの仕方を変え,動画の回復を試みる.動画を連続的に回復 することは計算コストが非常に高くなり,通常のフィルター処理を施すだけでも大変な作 業になってくる.従って,本研究においても現在のところは実時間においての回復ではな く,すでに撮影された動画をどれだけ高い精度で回復できるかどうかを確かめる事を目的 として実験を行った. 6.1 フレーム処理 動画の復元が困難な理由としては,フレームごとに処理系統が必要になってくるという 点である.具体的には,動画の 1 フレームを抜き出し,そのフレームに対して処理を施 し,次のフレームへ進むと言う作業が必要になってくる.また,もう 1 つの回復を困難に させている原因としては,画像のサイズである.フレームを抜き出して次のフレームへの 作業を連続的かつスムーズに行えたとしても,一枚一枚の画像サイズの大きさが画像回復 を難しくしている.本研究ではフレーム間の補正を行うことと,適応マスクの適用をする ことにより,セグメント化した小さな領域だけから劣化関数を推定できる事が出来るとい う利点から動画の回復を行う. 80 第 6 章 動画の復元 まず,フレーム間の補正であるが,こちらに関しては単純に前フレームの PSF の情報 を現フレームの初期予想 PSF として用いることにより,おおよその PSF の形がわかって いる状態から数値解析を進めることが可能であり,これにより第 2 フレームからの計算時 間をかなり減少することが可能となる.光学系に関して同じ位置にある物体からのインパ ルス応答はフレームの遷移という短い時間の間に大きく変化することはないと言う仮定の 下,ブラインド回復を行う.さらに,適応マスクを用いることにより画像全体の数値計算 を行う事無く,フレーム画像の一部を切り取ったものから PSF を求める事が可能である. 本研究では,適応マスクを用いることによるこの 2 つの利点を生かして動画のフレーム画 像の回復を行った. (a) (b) 図 6.1 動画の撮像系.2 通りの方法で動画を撮影する.(a) CCD カメラを物体に対 して垂直に移動させる. (b) CCD カメラに首振り運動をさせる. 6.2 実験 図 6.2 81 本論文において動画の回復に用いたセグメント.実線で囲まれた部分をセグメ ントとして用いる.劣化画像フレームに対して 5 カ所を固定して各々のセグメントに 対して PSF の推定を行う. 6.2 実験 6.2.1 実験方法 本研究では 2 通りの方法で動画を撮影し,各々のフレーム画像を回復することにより動 画の画像劣化の回復を行った.動画の撮像は図 6.1 に表されているとおりであり,1つ目 の方法としては,カメラを物体に対して垂直に動かす,2つ目の方法がカメラを首振りさ せた様な状態を想定した.撮影に用いた物体はテストチャートである.抜き出すセグメン トは図 6.2 に図示した. 82 第 6 章 動画の復元 実際に撮影した動画からいくつかフレーム画像を抜き出したものが図 6.3 と図 6.4 であ る.これらのフレーム画像から図 6.2 に従い5カ所のセグメントを抜き取り,抜き取られ た5カ所のセグメント情報から PSF の推定を行った. 6.2.2 PSF の推定 PSF の推定には各フレームにおいて遺伝アルゴリズムを組み込む.5 カ所のセグメン トからそれぞれ推定された PSF(h0 ∼h4 ) と得られた 5 つの推定 PSF の平均 (h5 ) を求め る.h0 ∼h5 を遺伝子として,GA に組み込む.ただし,突然変異確率は 0 とし,単純に 交叉・選択のみを行う.交叉方法は単純交叉,選択方法はエリート選択による. ・交叉:h0 ∼h5 からランダムに 2 つずつ選び単純交叉することにより h6 ∼h9 の 4 つの遺 伝子を新たに作成する. ・選択:各々の回復画像フレームと 10 個の PSFh0 ∼h9 をコンボリューションしたものと 劣化画像フレームの間のエネルギーを求め,各セグメントにおける PSF をエリート選択 により選ぶ. PSF の GA 過程は,マスクの更新処理と同時に行われる. 求められた PSF は次のフレームの初期予想 PSF として使用され,第 1 フレームでは解 析に時間がかかるものの,第 2 フレームからは第 1 フレームの大体 2∼3% ほどの時間で 解析が可能であった. 図 6.5 と図 6.6 にそれぞれ図 6.3 と図 6.4 を回復したフレームの画像を示す.また参考 として,図 6.3 の状態で10分後まで数値解析を行った際の復元像を図 6.7 に載せておく. 6.2 実験 図 6.3 83 C G D H E I F J 動画の劣化画像フレームを抜き出したもの.CCD カメラはテストチャートに 対して垂直に動かしている. (a) は初期フレームであり,(b) から (h) はそれぞれ (a) に続くフレームを 6 フレームごとに抜き出したものである.それぞれ 30FPS で撮影 された動画から 6 フレームごとに抜き出しているので,抜き出したフレームの間隔は 200ms である. 84 第 6 章 動画の復元 C G D H E I F J 図 6.4 CCD カメラを首振りさせた時の動画から劣化画像フレームを抜き出したもの. (a) から (h) のフレームは図 6.3 と同じ方法で抜き出されている. 6.2 実験 85 C G D H E I F J 図 6.5 図 6.3 の劣化画像フレームを回復したもの. 86 第 6 章 動画の復元 C G D H E I F J 図 6.6 図 6.4 の劣化画像フレームを回復したもの. 6.2 実験 図 6.7 10 分間画像を撮り続けた時の最終的な回復画像.フレーム数で言うと 18000 フレーム目の画像である.ただし,フレームの抜き出し間隔は 15 フレームごととして いる. 87 88 第 6 章 動画の復元 MTF frequency[1/pixel] 図 6.8 図 6.7 の回復像を得られた時の PSF の定量的な評価.規格化された MTF と 空間周波数のグラフ. 6.3 まとめ 89 6.3 まとめ 適応型マスクを用いることにより従来手法では困難であった動画の回復を行うことに成 功した.フレーム内ではセグメントを 5 カ所抜き出し,各セグメント内で PSF を推定す る過程はハイブリッド AD/SA 法を,セグメント間での PSF の推定には GA を用い PSF の最適化を行った. フレーム間の推定としては単純に前フレームの PSF を初期予想 PSF として代入する ことにより大幅な計算コストの削減に成功した. また,定量的な評価として x 方向の規格化された MTF を求めた.規格化された MTF は,以下の式の通り. ¯ ¯ ¯ F̃ (u, 0) + α ¯ ¯ ¯ M T Fn (u) = ¯ ¯. ¯ F (u, 0) + α ¯ (6.1) ただし,F̃ は最終的に求められた回復画像,F はベストピント画像,α は MTF の値が無 限大または 0 にならないようにするための微小な値を持ったパラメーターである.図 6.8 のグラフからもわかるとおり非常に高精度な回復像が得られていることがわかる. 91 第7章 結論 適応マスクという考え方を提案することにより,劣化画像のブラインド復元という分野 において非常に大きな貢献をした.適応マスクはブラインド・デコンボリューションの中 で用いられる手法ではあるが,その応用範囲としては,シフトバリアントな系における劣 化過程における劣化像や動画における劣化フレームの回復を行うことに成功した. 本研究における適応マスクアルゴリズムは,根底となるブラインド・デコンボリュー ションのアルゴリズムの種類にかかわらず単なるマスク処理として表すことが出来るの で,ブラインド・デコンボリューションの技術の成長と共に,その応用分野も広がってい くものと思われる.ブラインド・デコンボリューション法は,光学系にかかわらずロバ ストな回復を行うことが強みであるので,その応用分野も非常に広い範囲におよぶであ ろう. 全体的な流れとしては,第 1 章において,本研究の目的を示した. 第 2 章では,画像の PC 上での表現方法や基本的な画像処理について例を挙げた. 第 3 章では,従来行われているブラインド・デコンボリューション法の例を挙げ,その 問題点・解決すべき点を抑えて適応マスクの必要性を述べた.ゼロシートを用いることに より,コンボリューションによって得られた画像は唯一の物体と PSF に分離できること を説明した.また,ボルツマンマシンを用いたブラインド・デコンボリューション法のシ 92 第 7 章 結論 ミュレーションを行いブラインド・デコンボリューション法の実用性を示した. 第 4 章では,実際に光の移り込みが多い画像を準備し,計算機シミュレーションを行っ た.適応マスクは通常の窓処理やフィルター処理と違い,その時点で得られている推定 PSF を元にマスク処理を行うことにより,余分な情報を取り除く作用がある.また,その ときに得られている推定 PSF からマスクが作成されるので,ハニング窓やハミング窓な どのように縁の部分の情報を全く無いような状態にする様な窓処理ではないため,周波数 分解能はそれほど下がることなく,ROI 内部からの情報を多く消滅させてしまうことも 防いでいる.光の移り込みが多いと言うことは,適応マスクを適用しないことにはブライ ンド・デコンボリューション法を当てはめるのは困難であり,適応マスクを用いて PSF を正しい PSF に固定した上での SA 法と,適応マスクを適用し SA 法を用いた復元を比 較することにより適応マスクの有用性を実証した.また,コンボリューションによる劣化 過程 (シフトインバリアントな劣化過程) において,トランケートする部分を 1 箇所と 2 箇所の場合を比較した.セグメントの数を 1 つから 2 つにすることにより,計算時間はも ちろん,PSF の推定の精度に挙げることにも成功した.この点においてもセグメントか ら PSF を推定するために必要な適応マスクの効果が高く見られる.シミュレーションに 対し,実際に撮影された画像からの回復も試みた.撮影に用いられたのは,一眼レフデジ タルカメラであり,物体としては基本的なテストチャート画像を用いた.実験を行うに関 して,ノイズの影響を無視することは出来ないのでノイズ除去フィルターをかけるが,こ のとき用いるフィルターはデコンボリューションを妨げないようにするためにメディアン フィルターを用いている.実験結果からもわかるとおり,非常に高精度な画像復元を行う ことが可能であることが実証された. 第 5 章においては,第 4 章で紹介した適応マスクの応用としてシフトバリアントな PSF による劣化像の回復を行った.シフトバリアントな劣化過程の場合,通常はデコン ボリューションで行うことが不可能であるが,適応マスクを用いることによりデコンボ リューション技術を応用することが可能となり,計算の高速化,および精度の向上を期待 することが出来る.実際にシミュレーションにおいては非常に高精度な画像の復元を行う 93 ことに成功した.また,簡単な例として,前後離れた場所においた物体を同時に撮影した 劣化画像の回復を行った.前後に物体をずらすことにより,各々の物体に対する PSF を 変化させることによりシフトバリアントな PSF による劣化を表現している.実験として, 一眼レフデジタルカメラで撮影したもに,携帯電話に付属の CCD カメラで撮影したもの の復元を行った.一眼レフデジタルカメラで撮影された劣化像に関しては,非常に良い結 果が得られており,前後にずらされた物体の両方を高精度で復元することに成功した.対 して,携帯電話付属の CCD カメラを用いて撮影された劣化像に関しても良い結果が得ら れているが,携帯電話内で JPEG に変換されてしまっているため,JPEG 変換に依存す るブロックノイズが多くの領域で見られる.このため,CCD カメラモジュールそのもの の性能ではないところで画像復元に関するボトルネックが発生してしまっているため,簡 単に性能を測ることが困難になってしまっている.この研究に関してはカメラモジュール 自体を入手し劣化像を撮影し,その画像から回復する事が将来的な展望となってくるであ ろう. 第 6 章では,別の切り口からの適応マスクの応用例として動画を回復することを目的 として研究を行った.本論文で行われた実験はごく簡単な系に対してのものではあるが, 動画へ応用する基礎としては必要なものである.物体に対して CCD カメラを垂直方向お よび水平 (に近い) 方向へ連続的に動かすことにより,ピントぼけやモーションブラーを シミュレートした.どちらの劣化動画に対してもフレームごとにブラインド・デコンボ リューション法を適用している.ただし,各フレームからセグメントを抽出してブライン ド・デコンボリューションを適用することは非常に計算的なコストが高くなり,そのまま の状態で計算を行うことが困難になってくる.そのため,一つ前のフレームで得られた推 定 PSF を用いることにより計算時間の短縮を行っている.また,第 4 章で実証された複 数のセグメントによる高精度な回復を元に,セグメントの数を 5 個にして計算を行ってい る.5 個のセグメントに対して,各々 PSF を推定した後,マスクの更新のタイミングで 遺伝アルゴリズムを用いることによりさらなる高性能化を目指した.遺伝アルゴリズム過 程では,交差することにより 10 個まで遺伝子を増やし,その中からエリート選択をする 94 第 7 章 結論 ことによってよりよい PSF だけを次の世代へつなげるようにした.この過程を用いるこ とにより SA 法だけに頼るデコンボリューションよりも高速に計算を行うことを可能とし ている. 例えば,天体画像から始まって,顕微鏡で見るミクロな世界,また医療分野においても 簡単な光学系で組むことの出来るカプセル型内視鏡などが挙げられるし,身近な分野では デジタルカメラにおける画像処理系への応用,今より簡単な光学系による携帯電話付属の カメラなどである. 適応マスクの考え方は,デコンボリューションにおいて非常に有用なものであり,デコ ンボリューション技術はこれからも発展していく物であるので,それに伴って必要性も上 がってくる技術になる可能性を秘めている技術である. 本研究においては非常に簡単な作成方法である適応マスクを用いることによってブライ ンド・デコンボリューション法に新たな可能性を生み出した. 95 謝辞 本論文を作成するにあたり,学部在学中より長い間終始ご指導いただき,研究を行う場 と機会を与えて下さった早稲田大学先進理工学部応用物理学科小松進一教授に心より感謝 の意を表します.共同合宿などにおいて適切なご指導とご助言を頂いた早稲田大学先進理 工学部応用物理学科大頭仁教授,鵜飼一彦教授,また,論文作成に際してご指導とご助言 を頂いた早稲田大学先進理工学部応用物理学科橋本周司教授,森島重生教授に深く感謝の 意を表します. 本研究は日本学術振興会特別研究員として 2003 年から 2005 年の間,また早稲田大学 理工学部応用物理学科助手として 2005 年から 2008 年の間に行った研究であり,両研究 機関には研究を行う機会を頂いたことに深く感謝の意を表します.同時に両機関からの助 成により本研究がなされた事を併せてここに記します. また,経済産業省特許庁第四部画像処理部屋の皆様には画像処理に関する基本的な事か ら視野を広げた分野に関してご指導頂いたことに深く感謝の意を表します. 早稲田大学先進理工学部小松進一研究室,大頭仁研究室,鵜飼一彦研究室の皆様からは 学ぶことも励まされることも多く非常にお世話になりました.特に服部雅之研究員,日本 女子大学理学部小川賀代講師の両氏には学部時代から研究や進学,仕事についての事まで 様々な相談にのって頂いたことに心より感謝の意を表します. 最後に,9 年と言う長い学生生活にも関わらず,常に協力下さり,支えになって頂いた 父 規生,母 史江,姉 充子,由希子に心より感謝の意を表します. 97 参考文献 [1] J. W. Goodman, Introduction to Fourier Optics, 3rd ed., Roberts & Company Publishers (2004/12). [2]「光とフーリエ変換」谷田貝豊彦,朝倉書店 (1992/5). [3] L. B. Lucy, ”An iterative technique for the rectification of observed distributions,” Astron. J. 79, 745 (1974). [4] W. H. Richardson, ” Baysian-based iterative method of image restoration,” J. Opt. Soc. Am. 62, 55 (1972). [5]「科学計測のためのデータ処理入門 第 2 版」南茂夫 監修,河田聡 編著,CQ 出版社 (2002/5). [6]「ディジタル画像処理」CG-ARTS 協会 (2006/3). [7]「C 言語で学ぶ実践画像処理」井上誠喜 他,オーム社 (1999/12). [8] R. G. Lane and R. H. T. Bates, ”Automatic multidimensional deconvolution,” J. Opt. Soc. Am. A4, 180 (1987). [9] R. Fletcher and C. M. Reeves, ”Function minimization by conjugate gradients,” Comp. J. 7, 148 (1964). [10]「共役勾配法」戸川隼人,教育出版 (1977/7). [11] R. A. Fisher, ”On the mathematical foundations of theoretical statistics,” Philosophical Transactions of the Royal Society, A, 222, 309 (1922). [12]「遺伝アルゴリズムの理論」John H. Holland (嘉数侑昇 監訳),森北出版 (2005/11). 98 参考文献 [13]「ニューロ・ファジィ・遺伝的アルゴリズム」萩原将文,産業図書 (1994/7). [14]「ニューラルネットワークの基礎と応用」馬場則夫 他,共立出版 (1994/9). [15] T. Kohonen, ”The Neural Phonetic Typewirter,” IEEE computer, 21, 11 (1988). [16] T. Kohonen, G. Barna and R. Chrisley, ”Statistical pattern recognition with neural networks: benchmarking studies,” Proc. ICNN, I, 61 (1988). [17] F. Rosenblatt, ”The Perceptron: A Probabilistic Model for Information Storage and Organization in the Brain,” Psychological Review, v65, 386 (1958). [18] J. J. Hopfield, ”Neural network and physical systems with emergent collective computational abilities,” Proc. National Academy of Sciences of the U.S.A. 79, 2554-8 (1982). [19]「乱数とモンテカルロ法」宮武修,脇本和昌,森北出版 (1978/8). [20] D. H. Ackley, G. E. Hinton and T. J. Sejnowski, ”A learning algorithm for Boltzmann machines,” Cognitive Science, 9, 147 (1985). [21] S. Kirkpatrick, C. D. Gelatt and M. P. Vecchi, ”Optimization by Simulated Annealing,” Science 220, 671 (1983). [22] J. R. Fienup, ”Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758 (1982). [23] G. R. Ayers and J. C. Dainty ”Interative blind deconvolution method and its applications,” Opt. Lett. 13, 547 (1988). [24] S. Komatsu and J. C. Dainty, ”Blind deconvolution of a blurred picture using an iterative algorithm,” Appl. Opt. Dig., 295 (1990). [25] N. Suetake, M. Sakano and E. Uchino, ”Gradient-Based Blind Deconvolution with Phase Spectral Constraints,” Opt. Review 13, 417 (2006). [26] D. Ghiglia, L. Romero, and G. Mastin, ”Systematic approach to two-dimensional blind deconvolution by zero-sheet separation,” J. Opt. Soc. Am. A10, 1024 (1993). 99 [27]「工学基礎 ラプラス変換と z 変換」原島博,堀洋一,数理工学社 (2004/10) . [28] R. H. T. Bates, B. Quek, and C. Parker, ”Some implications of zero sheets for blind deconvolution and phase retrieval,” J. Opt. Soc. Am. A 7, 468 (1990). [29] B. C. McCallum, ”Blind deconvolution by simulated annealing,” Opt. Commun. 75(2), 101 (1990). [30] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller, ”Equations of State Calculations by Fast Computing Machines,” J. Chem. Phys., 21, 1087 (1953). [31] T. J. Holmes, ”Blind deconvolution of quantum-limited incoherent imagery: maximum-likelihood approach,” J. Opt. Soc. Am. A 9, 1052 (1992). [32] T. J. Schulz, ”Multiframe blind deconvolution of astronomical images,” J. Opt. Soc. Am. A 10, 1064 (1993). [33] E. Y. Lam and J. W. Goodman, ”Iterative statistical approach to blind image deconvolution,” JOSA A, 17, 1177 (2000) . [34] I. Numata and S. Komatsu, ”Improved Blind Deconvolution Algorithm Using a Sub-Band Filter Bank,” Technical Report (Advance Research Institute for Science and Engineering, Waseda Univ.), No. 2005-4 (2001). [35] R. J. Steriti and M. A. Fiddy, ”Blind deconvolution of images by use of neural networks,” Opt. Lett. 19, 575 (1994). [36] H. Tanaka and S. Komatsu, ”Blind deconvolution algorithm by use of quasi binary Hopfield neural network,” Technical Report (Advance Research Institute for Science and Engineering, Waseda Univ.), No. 2005-5 (2001). [37] 石原信人,小松進一,“ニューラルネットワークモデルを用いた劣化像のブラインド 回復法,”O plus E,28 巻 3 号,258 (2006). [38]「遺伝アルゴリズムを用いたブラインドデコンボリューション」石原信人,早稲田大 学理工学部応用物理学科 卒業論文 (2000). 100 参考文献 [39] Y. Chen, Z. Nakao, K. Arakaki and S. Tamura, ”Blind Deconvolution Based on Genetic Algorithms,” IEICE Trans. Fund. 80-A, 2603 (1997). [40] T. Enokura, Y. Chen, Z. Nakao and S. Tamura: Proc. IEEE Int. Conf. Systems, Man and Cybernetics, 127 (1999). [41] 山口昌一郎,石原信人,小松進一,”免疫アルゴリズムを用いたブラインド・デコン ボリューション,” Proc. Optics & Photonics Japan 2006,392 (2006). [42] J. N. Caron, N. M. Namazi and C. J. Rollins, ”Noniterative blind data restoration by use of an extracted filter function,” Appl. Opt. 41, 6884 (2002). [43] N. Wang, Y. Chen, Z. Nakao and S. Tamura, Parallel-Distributed Blind Deconvolution Based on a Self-Organizing Neural Network,” Appl. Opt. 38, 4345 (1999). [44] N. Miura, K. Ohsawa and N. Baba, ”Single-frame blind deconvolution by means of frame segmentation,” Optics Letters, 19, 695 (1994). [45] N. Miura, ”Blind deconvolution under band limitation,” Opt. Lett. 28, 2312 (2003). [46]「ディジタル信号処理」 辻井重男,鎌田一雄,昭晃堂 (1990/04). [47] N. Ishihara and S. Komatsu, ”Dynaic Masking of Truncated Image Segments for Blind Recovery of Blurred Images Based on Simulated Annealing,” Opt. Eng., 43, 1177 (2004). [48] 小出温子,小松進一, ”シフトバリアントな psf によって劣化した画像の復元−行列 演算反復法−,” 光学 26, 280 (1997). [49] Nobuhito Ishihara and Shinichi Komatsu, ”Blind Recovery for Degraded Image Blurred with Shift Variant PSF: Series Segmenting Method,” Technical Digest ICO’04 (Tokyo, Japan, Jul. 2004), 445 (2004) [50] 石原信人,小松進一,”連続分割法を用いたシフトバリアント劣化像のブラインド回 復 −実験的検証−,” Optics Japan 2004 (日本光学会年次学術講演会) 予稿集, 160 101 (2004) 研究業績 学術論文 • Nobuhito Ishihara and Shinichi Komatsu, ”Dynamic masking of image deconvolution for moving pictures,” Jpn. J. Appl. Phys., Vol. 46, Part 1, No. 12, 7763-7767 (2007) • Nobuhito Ishihara and Shinichi Komatsu, ”Dynaic Masking of Truncated Image Segments for Blind Recovery of Blurred Images Based on Simulated Annealing,” Opt. Eng., Vol. 43, No. 5, 1177-1182 (2004) 国際会議論文 • Nobuhito Ishihara and Shinichi Komatsu, ”A novel method of blind de-blurring for consecutive frames of a motion picture,” Technical Digest FiO’07 (San Jose, USA, Sep. 2007), JSuA21.pdf (2007) • Nobuhito Ishihara and Shinichi Komatsu, ”Dynamic Masking of Image Deconvolution for Micro Order Moving Pictures,” Technical Digest MOC’06 (Seoul, Korea, Sep. 2006), 138-139 (2006) • Nobuhito Ishihara and Shinichi Komatsu, ”Blind recovery of Blurred Image Using Adaptive Masking Method: Experimental Results,” Proc. MOC’04 (10th Microoptics Conference, Jena, Germany, Sep. 2004), L 59.pdf (2004) • Nobuhito Ishihara and Shinichi Komatsu, ”Blind Recovery for Degraded Image Blurred with Shift Variant PSF: Series Segmenting Method,” Technical Digest ICO’04 (Tokyo, Japan, Jul. 2004), 445-446 (2004) • Nobuhito Ishihara and Shinichi Komatsu, ”Blind Recovery of Truncated Blurred Image Using Adaptive Masking Method,” Proc. WISICT 2004 (Cancun, Mexico, Jan. 2004), 268-272 (2004) 総説 • 石原信人,小松進一,“ニューラルネットワークモデルを用いた劣化像のブライン ド回復法,” O plus E,28 巻 3 号,258-262 (2006) 講演 • 石原信人,小松進一,“動画のための連続分割画像回復法,” Optics & Photonics Japan 2006 (日本光学会年次学術講演会・日本分光学会秋期講演会) 予稿集, 394-395 (2006) • 山口昌一郎,石原信人,小松進一,“免疫アルゴリズムを用いたブラインド・デコ ンボリューション,” Optics & Photonics Japan 2006 (日本光学会年次学術講演 会・日本分光学会秋期講演会) 予稿集,392-393 (2006) • 石原信人,小松進一,“ニューラルネットワークモデルを用いた連続分割法による シフトバリアント劣化像のブラインド回復,” Optics Japan 2005 (日本光学会年 次学術講演会) 予稿集,136-137 (2005) • 石原信人,小松進一,“連続分割法を用いたシフトバリアント劣化像のブラインド 回復 −実験的検証−,” Optics Japan 2004 (日本光学会年次学術講演会) 予稿集, 160-161 (2004) • 松澤重信,石原信人,小松進一, “ブラインド・デコンボリューションのための自己 組織的アルゴリズム,” Optics Japan 2004 (日本光学会年次学術講演会) 予稿集, 162-163 (2004) • 石原信人,小松進一,“適応型マスク法を用いたブラインドデコンボリューション による像回復, ” 第 64 回応用物理学会学術講演会講演予稿集, No. 3, 905 (2003) • 松浦大和,石原信人,田島由理,小松進一,“SA 法に基づくシフトバリアント劣化 画像のブラインド回復, ”第 64 回応用物理学会学術講演会講演予稿集, No. 3, 905 (2003) • 石原信人,小松進一,“遺伝アルゴリズムを用いたブラインドデコンボリューショ ン,” 第 61 回応用物理学会学術講演会講演予稿集, No. 3, 880 (2000)