画像復元における正則化パラメータの推定 利用統計を見る

論 文

画像復元における正則化パラメータの推定

渡部知樹 大木真 周欣欣 橋口住久

（平成8年8月31日受理）

Estimation of the Regularization Parameter

for Image Restoration

SatokiWATANABE MakotoOHKI XinxinZHOU SumihisaHASHIGUCHI

Abstract   In the image restoration using the regularization method， it is important to choose the reg． ularization parameter apPropriately・However， the regularization parameter is usuaUy unknown apriori． Although there are sever訓estimation methods for the regularization parameter such as cut and try method and the statistical estimation method，they require a huge amount of computation．   This paper proposes an estimation methods for the regularization parameter using simple calculation． In order to reduce the computation， this method exploits the intermediate output of theCombined DFT／LDE（Discrete Fourier Transform／Linear Difference Equation）methods which is used to implement the image restoration by the 2．D recursive filter． Although the estimation accuracy of the ploposed method is not high， the quality of the restored image is near to that of the optimal one．

1 まえがき

 画像復元の目的は，ぼけやぶれなどにより 劣化した画像から，計算機処理などによって， 原画像にできるだけ近い画像を求めることで ある．再帰形フィルタを用いた画像復元法の 一つとして・正則化フィルタ（RegularizatiQn filter）1）を用いた復元法があり，その実現法と して，Combined DFT／LDE（Discrete Fourier Transform／Linear Difference Equation）法2） を用いて，2次元のスペクトル因数分解を1次 元のスペクトル因数分解に帰着する方法が提案 されている3）．  正則化フィルタを用いて実際に画像復元を行 なうには，劣化画像から適切な正則化パラメー タを推定する必要がある．本研究では，Com− bined DFT／LDE法の中間結果として得られる 周波数スペクトルを利用して，大まかな近似と 簡単な計算によって正則化パラメータを推定す る方法を提案する．

2 正則化と正則化フィルタ

画像の劣化過程は以下の式でモデル化される． ＊電子情報工学科，Department of Electrical Engineering and Computer S cience 9（k・，k2）＝碓・，k・）・∫（k・，k、）＋n（le、，k2）       （1）

ここで，h（le1 ， k2）は画像の劣化を表わすイン パルス応答，∫（克1，k2）は原画像，g（た1，k2）は 観測画像，n（k1，k2）はノイズを表わす．ぼけの 周波数特性H（ω1，ω2）は2次元であるが・説 明の便宜上，その1方向からの断面のみを示 すと図一1のようになる．図からわかるように， H（ω1，ω2）は通常，零点を含む低域通過特性と なる．このH（ω1，ω2）の逆特性は図一2のよう になり，そのままでは安定な伝達関数とならな い．また，高域通過特性であるために，高周波 域のノイズが過剰に強調されてしまう．  そのため，図一3のように，逆フィルタの発 散しているところを安定になるように抑え，ま た，高域の利得を下げることによって，高域の ノイズの過剰な強調を抑える．これを正則化 という．正則化の程度は正則化パラメータλに よって調節される．図4は正則化パラメータを 変えて行った時の復元誤差のモデルである．正 則化パラメータが小さすぎる場合，復元フィル IHI 0      ω     図一一1 Hの周波数特性 Fig．1 Frequency Characteristics of H． IH−11  0       ω     図一2 H−1の周波数特性 Fig．2 Frequency Characteristics of H−−1． 1助   0      ω  図一3 正則化フィルタWの周波数特性 Fig．3 Frequency Characteristics of Regular− ization］Filter W． Error 0 λ 図一4 正則化パラメ・・…一タに対する復元誤差 Fig．4 Restoration Error v．s． for Regulariza− tion Parameter． タの高域利得が高くなっているので、ノイズが 強調されてしまい，ノイズ拡大誤差が大きくな り，これによって，復元誤差も大きくなってし まう．また，このノイズ拡大誤差は画像全体に 作用する．逆に，正則化パラメータが大きすぎ る場合，画像に滑らかさを強調するあまりエッ ジ付近にぼけが生じるために，あまり画像が復 元されず，正則化誤差が大きくなるので，やは り復元誤差が大きくなる．  このように，ノイズ拡大誤差と正則化誤差は トレードオフの関係となり，正則化パラメータ には最適値が存在することになる．劣化画像か ら適切な正則化パラメータを求める方法として は，統計的に推定する方法4）や試行錯誤による 方法が提案されているが，多大な計算量が必要 となる．これに対し，本研究では，大まかな近 似と簡単な計算によって，自動的に正則化パラ メータを推定する方法を提案する．

平成8年12月  式（1）で表される観測画像g（kl，k2）から，原 画像∫（k1，k2）を推定するために，式（2）の評価 関数を最小にするように画像復元フィルタを設 計する．ここで，∫（kl，k2）は復元画像である．

    

ξ・＝Σ［h（k・， k2）・∫（ゐ・，た・）−9（ゐ・，先・）］2   た1，k2＝一’◎o

＋λ￡卜⑭・∫⑭］2（2）

   先11た2＝一◎o 第1項は原画像と復元画像の誤差の量を表わ し，第2項は復元画像に含まれるノイズ成分の 量を表わしている．λは正則化パラメータであ り，また，c（先1，克2）は正則化演算子である．  正則化演算子としては，ラプラシアンフィル タなどの高域通過フィルタを用いるのが一般的 である．これは，以下のような理由による． 画像の電力スペクトルは低周波域に集中して おり，また画像の劣化特性も多くの場合低域通 過特性を持っので，劣化画像の高周波域成分は ノイズ電力がほとんどであると考えられる．し たがって，復元画像におけるノイズ成分の量を 評価するためには，復元画像の高周波成分を取 り出してやればよい．  式（2）の偏微分を0とおいて∫について解き， フーリエ変換すると，画像復元フィルタの周波 数特性は， なければならないので，計算量，計算時間とも に膨大になってしまう．また，正則化パラメー タの推定の研究もいくっか行なわれている4） が，これらは，反復的な行列計算による復元法 に統計的な推定を組み合わせたものであり，や はり，多大な計算量を必要とする．ここでは， 大まかな近似と簡単な計算，また，Combined DFT／LDE法から自動的に出てくる周波数ス ペクトルを用いて，正則化パラメータの推定を 行なう．そこで，本研究では，古典的なウィー ナーフィルタを基準として考え，ウィーナーフィ ルタに類似した特性が得られるように，劣化画 像から自動的に正則化パラメータを推定する方 法を示す．  まず，ウィーナーフィルタの周波数特性を次 式で表す。 w（ω1，ω2） H＊（ω・，ω2） 1∬（ω・，ω2）12＋5九（ω・，ω2）／5∫（ω、，ω2）        （4） Sf（ω1，ω2）およびSn（ω1，ω2）は原画像および ノイズの電力スペクトルである．また，正則化 フィルタの周波数特性は， 171zλ（ω、，ω2）＝ 王r＊（ω1，ω2） IH（ω・，ω2）12＋λ10（ω・，ω2）12       （5） vvλ（ω1μ2）＝ H＊（ω・，ω2） となる．    1∬（ω、，ω2）12＋λ10（ω1μ2）12        （3） このとき，＊は複素共役を表わす．

3 正則化パラメータの推定

3．1 正則化パラメe−m一タ推定式の原理  試行錯誤によって正則化パラメータを求める 場合，まず，適当に任意の正則化パラメータの 値を正則化フィルタのλに入れ，画像を復元し， 原画像と誤差を比較するということを，誤差が 小さくなるように反復的に何度も復元しなおさ であるから，式（4）と式（5）が等しくなるため には，正則化パラメータが x＝ 5・（ω・，ω・）   10（ω・，ω2）125∫（ω・μ2） （6） となれば良い．式（6）の分母は原画像の電力ス ペクトルを含んでいるが，実際には，原画像の 情報は得られないので原画像の電力スペクトル は使えない．そこで，式（1）を電力スペクトル で表すと， 5、（ω・，ω・）＝1∬（ω・，ω・）125∫（ω・，ω、）＋Sn（ω、，ω、）        （7）

であるから，これを変形して， 5∫（ω・μ2）

   1

「∬（ω、，ω、）1・［Sg（ω・，ω・）−Sn（ω・，ω・）】 （8） とする．ここで，ノイズの電力スペクトルは， 劣化画像の電力スペクトルよりも十分小さいか ら，Sg（ω1，ω2）一一 Sn（ω1，ω2）2t Sg（ω1，ω2）と近 似すると，         1 Sf（ω・，ω・）＝ hH（ω、，ω、）1・Sg（ω・，ω・）（9） これを式（6）に代入してSf（ω1，ω2）を消去する と正則化パラメータの推定式は， x＝ Sn（ω・，ω・）   1麗i：li；：IPsg（ω・，ω・） （10） 式（10）の分母のHω1μ、）は・図一2のように 発散する恐れがあるので，そのまま用いること ができない・そのため・おおよそ珂㍍｝こ近 い形で安定であることが保障されている特性 H。h（ω1，ω2）を作り，これで置き換える．ここ で・H（wl；it；2）・H・h（ω・lw・）｝ま，2次元の周波数 特性である．  本研究では，Combined DFT／LDE法を用 いて画像復元を行なっているが，Combined DFT／LDE法では，克2方向にのみ離散フーリエ 変換した周波数スペクトルSg（kl 1ω2）が中間結 果として得られる．そこで、2次元の周波数ス ペクトルを直接計算する代わりに、この周波数 スペクトルSg（kl， w2）を用いることを考える． ピンボケなどによる劣化は円対称であるので， 1次元の周波数スペクトルを用いても十分な推 定が可能である．手ぶれなどによる非対称な劣 化にっいては今後の課題である．劣化画像のス ペクトルとしてSg（ゐ1，ω2）を用いるとすると， 0・Hah等も1次元で考えれば良い．式（10）を 近似計算すると，以下のようになる． x＝  s九（k1，ω2）       （11）   1∬。九（ω2）1210（克・，ω2）1259（克・，ω2） 式（11）のように推定式は導出される．  また，安定な近似フィルタHah（ω2）の具体 的な形は次の通りである． H。h（ω、）＝1＋α［1−（1／2＋1／2…（ω、））N／2］        （12） 式（12）のdZV’は劣化特性H（ω1，ω2）の次数，α はHahの高周波域における振舞を調節するパラ メータである．シミュレーションの結果から， 一様なぼけの特性の場合，αはぼけのインパル ス応答の振幅の逆数の平方根で表され，特に， 円対称のぼけで次数が低次である場合には， α＝   N2／2十N十1 （13） とすることにより，適切なHah（ω2）が得られ ることが示された．

3．2 ノイズの電力スペクトル

 ノイズの電力スペクトルは，劣化画像から直 接得ることができない．そこで，まず，ノイズ は白色ノイズであると仮定する．これによって， ノイズは周波数に関係なく一定であるとみなせ る．ぼけの周波数特性は，低域通過特性である ので，図一5のように，原画像にぼけのみをかけ た劣化画像の電力スペクトルは，原画像の電力 H 日 竃 合 8 き

A

0 frequency _T    図一5 電力スペクトル特性 Fig・5 Power Spectrum Characteristics。

平成8年12月 第47号 スペクトルSfと比較して，高域周波数成分が ほとんどなくなり，Sbのようになる．そこに， ノイズSnが加わり，劣化画像の電力スペクト ルSgとなる．図一・5のSgのように，劣化画像の 電力スペクトルの高周波数成分は，画像の電力 スペクトルはほとんどなく，ノイズの電力スペ クトルだけであると考えられる．また，劣化画 像の電力スペクトルの高周波数成分は，ほぼ一 定であるので，ノイズの電力スペクトルは，劣 化画像の電力スペクトルの高周波域におけるレ ベルSghとして近似できる．

3．3 正則化パラメータの推定式

以上から，式（11）は，

入＝     Sgh

      （14）   IH。九（ω2）121C（k・，ω2）12Sg（k・，ω2） となる・式（14）のSghは劣化画像の電力スペ クトルの高域成分，Sg（kl1ω2）は彪方向のみに 離散フーリエ変換した観測画像のスペクトル， 0（ω1，ω2）は正則化演算子を表す．本研究では， 式（14）を用いて正則化パラメータの推定を行 なう．  具体的には，次のような手順で推定を行う． 1．劣化画像を横方向に離散フーリエ変換した  画像Sg（k1，ω2）から適当に1行分のデー  タを選び，高域成分（仮定によりω2に対  して平坦である）の平均をとって求めた定  数値を式（14）の分子Sghとする． 2．同じ1行分のデータに対してIHah（ω2）12，  10（k1，ω2）12というフィルタをかける． 3．フィルタ結果はω2に対してほほ平坦とな  る（そうなるようにHah（ω2）を定めた）の  で，ω2方向に平均をとって求めた定数値  （仮にdと置く）を式（14）の分母とする． 4．叉＝竿によって正則化パラメータの推  定値を求める． 5．1．∼4．を画像の各行にっいて行ない，そ  の平均を正則化パラメータとして用いる．

4 シミュレーション

4．1 シミュレーション方法

 使用する画像は，大きさ256×256pixels，8 bit／pixelの白黒標準画像aerial， lenna， man− drillを用いてシミュレーションを行なった．ぼ けの特性hとして以下に示すような，ぼけの 大きさ（COC＝Circle Of Confusion）が7の一 様なぼけを用いた．また，Cは正則化演算子で ある高域通過形フィルタである．使用した復元 フィルタは・6次のものを2次の縦続形で実現 している． ん（k1，k2）＝ 0 0 0 0．04 0 0 0 c（kl，h2）＝

［÷

0    0   0．04   0    0 0   0．04  0．04  0．04   0 0．04  0．04  0．04  0．04  0．04 0．04  0．04  0．04  0．04  0．04 0．04  0．04  0．04  0．04  0．04 0   0．04  0．04  0．04   0 0    0   0．04   0    0

li：：剥

0 0 0 0．04 0 0 0 （15） （16）  観測画像には，SNR＝28dBとなるように白 色ガウスノイズを加える．ただし，観測画像の SNRは次式で定義される． SNR＝・・1・gva「1

f。蒜1豊豊芸age

4．2 画像復元結果

（17）  本研究による推定結果を表一1に示す．また， 本推定法と比較するために，MSEが最小に なるように試行錯誤で求めた正則化パラメー タλMMSEと，その復元誤差MMSE（Minimum MSE）を用いる．

ここで，MSE（平均2乗誤差）の定義は，     M芸1嘗1（∫（k1，k2）一∫（le1， le2））2     k・＝妾 k2＝f

MSE＝

        （M−」）×（N−J）        （18） であり，∫（k1，k2）は原画像，∫（le1 ， k2）は復元 画像，M×Nは画像の大きさである．ただし， 過渡応答によって，画像の端の部分が正しく復 元されないので，この過渡応答の影響を除くた め，回り9画素づっ無視したものをMSEとし て用いている．ここで，Jは，フィルタの次数 を示す． 表一一1λの推定値とMSEが最小になるλの比較 Table−1  Comparison between the E8ti． mated Value and the Optimal One． 画像名 λ

M51ヲ λMMSE

M5E

aerial _{0，033 108．7  0．019 106．3} lenna 0，038 79．9   0．13 71．8 mandrill 0，027 199．6  0．012 190．2  シミュレーションの例として，mandrillの原 画像，劣化画像，MSEが最小の復元画像，λの 推定による復元画像をそれぞれ図一6，図一7，図 一8，図一9に示す．

4．3 シミュレーション結果の考察

 本推定法によって求められた正則化パラメー タλと，MSEが最小になるように試行錯誤で 求めた正則化パラメータλMMSEを比較すると， λの推定精度は高くはなく，また，画像による バラツキも大きいことがわかる．しかし，MSE を比較するとあまり差がなく，復元画像を比較 しても，ほとんど同等の画像が得られている． これは，図4を見ればわかるように，正則化パ ラメータが最適値付近である場合，復元誤差の 特性はなだらかになっており，多少の推定誤差 はあまり影響してこないためである。  次にパラメータ推定の計算量を考察する．以

下では画像のサイズをM×Nとする．本推

定法の主な計算は，3．3の手順に示すように画 像の1行分に対する周波数域でのフィルタリン グであり，それを各行について行なっているの で，乗算回数のオーダーはMNの程度である．  他方，統計的推定法の場合4）は，近似的な推 定値を得るために，画素値を辞書式に並べたベ クトル（サイズMN）に対して行列の乗算を 行うので，そのままでは（MN）2のオーダーの 乗算回数が必要となる．何らからの高速算法を 用いたとしてもMN Iog2 MNのオーダーであ る．さらに，推定値が収束するまで推定を反復 的に繰り返す必要があり，反復回数が何回とな るかは画像に依存しておりあらかじめ知ること はできない．ただし，推定値の精度は本推定法 よりも優れている．  また，試行錯誤による方法は正則化パラメー タの値を適当に変えながら画像復元を繰り返す ものであり，パラメータ推定のために画像復元 そのものの計算量の何倍もの計算量を必要とす ることになるので，計算量の点でのメリットは ほとんどない．ただし，簡便であるので計算時 間に制限がないときには有用である．

5 まとめ

 本研究では・Combined DFT／LDE法によ る画像復元において，大まかな近似と簡単な計 算によって正則化パラメータを求める方法にっ いて述べた．  復元結果については，，円対称なぼけの特性の 場合，本方法によって推定した正則化パラメー タによる復元画像は，最適な正則化パラメータ による復元画像と比較して，視覚的にはほぼ同 等であり，また，MSEにも大きな差が認めら れないような結果が得られた．

参考文献

1）M．Unser， A． Aldroubi and M． Eden  ：‘‘Recursive Regularization Filters：  Design，Properties，and Applica一

平成8年12月 第47号   tion，，， IEEE Transactions on Pattern   Analysis and Machine Intelligence，   Vo1．13， No．3， pp．272−−277， March．   1991． 2）L．T． Bruton and N． R． B artley   ：“Application of Complex Filters   to Re aJize Three−Dimensional Com−   1）ined DFT／LDE Transfbr Func−   tions，，， IEEE Tra皿sactions on Cir−   cuits and Systems， Vo1．39， No．6，   pp．391−394， June．1992． 3）大隅升男，大木真，橋口住久：“Com−   bined DFT／LDE法による再帰形画像   復元フィルタの実現”，1994年電子情報   通信学会春季大会，A−118， p。189，1994   年3月。 4）Stanley J． Reeves：Optimal Space−   Varying Regularization in Iterative   Image Restoration  IEEE Trans−   actions on Image Processing，Vo1．3，   No．3， pp319−324， May．1994．    図一6 原画像（mandrill） Fig．6 0riginal Image（mandrill）． 図一8 MSEが最小の復元画像（λ＝0．012） Fig．8 Restored Image with Minimim MSE （λ：＝0．012）．     図一7 観測画像（SNR＝28dB） Fig．7 Blurred Image with Gaussian Noise （SNR＝28dB）． 図一9λの推定値を用いた復元画像（λ＝o．027） Fig．9 Restored Image Using the Estimated λ（λ＝0．027）．