論 文
画像復元におけるぼけのインパルス応答の推定
阿部雅文 大木真 周欣欣 橋口住久
（平成9年8月29日受理）
Estimation of the Impulse Response of Blurs
for Image Restoration
MasahumiABE MakotoOHKI XinxinZHOU SumihisaHASHIGUCHI
Abstract This paper proposes a estimation method for the impulseごesl）onse of the image degradation model from the degraded image． U皿der assumptio皿s that the noise is zero−mean white Gaussian noise and the degradation is the uniform horizontal blur， the proposed method can estimate the impulse response by counting the number of zero−crossing points of the power spectrum of the degraded image．
1 まえがき
画像として得られる情報には，カメラや顕微 鏡などの光学的レンズによって得るものや，テ レビジョンシステムなどによるものがある．と ころが，その画像が原場面をすべて正確に表し ているとは言い難い．例えば，写真を撮影する 場合，手ぶれなどにより画像がぼけてしまい，更 にその画像に雑音が加わり，原場面とは異なっ てしまう．このような画像は劣化画像と呼ばれ る．劣化画像は見苦しいだけでなく，対称物の 形状の把握や対称物に含まれる特徴や情報の抽 出を困難にする．そこで，観測された劣化画像 から計算機処理によって，元の画像にできるだ け近づけて必要な情報を読みとれるようにする 必要がある．このような処理は画像復元と呼ば れる， 画像復元には，復元フィルタを用いる．復元 フィルタを設計する際には，劣化に関わるぼけ の種類，大きさ，雑音の統計的性質（これらを劣 化モデルという）が既知でなくてはならない． 従って，劣化画像から劣化モデルを推定すると いう操作が別個に必要である．この操作は画像 同定と呼ばれる． ＊電子情報工学科，Department of Electrical Engineering and Computer Science 画像同定の手法としては，画像の生成モデル と雑音の統計的性質に適当な仮定を置いて最尤 推定などによって劣化モデルを推定する方法が 提案されている1）．しかし，この方法は非線形 の最適化問題を解くために多大な計算量を必要 とするので，一般民生用機器に用いるのは困難 である．一方，推定精度は劣るが簡便な推定法 としては電力スペクトルの特徴に注目する方法 がある．例えば，文献2）では，雑音は白色ガウ ス雑音，劣化は一様水平ぼけという仮定の下で 劣化画像の電力スペクトルの特徴から雑音分散 の推定を行なっている．ただし，この論文では， ぼけの大きさは既知と仮定している． 本研究では，文献2）と同じ仮定の下で劣化画 像の電力スペクトルの特徴からぼけの大きさを 推定する方法を提案する．2 画像復元
2．1 画像の劣化と復元
図一1のように，原画像fをイメージングシ1 テムHを通して観測すると，手ぶれや焦点ぼ けによりぼけを生じ，更に雑音nが加わること により，観測画像gは劣化する．図一2のように， 劣化した観測画像を復元フィルタWを通すな どの計算機処理によって，劣化した画像からぼf 9 原画像
H
劣化画像 イメージング システムn
雑音 図一1 画像の劣化モデル Fig．1 Degraded image model9
劣化画像 復元フィルタ 図一2 画像の復元モデル Fig．2 Restorated image model けや雑音を取り除き復元画像fを得ることが画 像復元の目的である．ここで復元フィルタを設 計する際に劣化モデルが既知である必要がある 1）3）．2．2 画像の同定問題
2．2．1 画像の同定問題 2．1より画像復元を行なうには，劣化に関わ るぼけの種類，大きさ，雑音の統計的性質（劣化 モデル）が既知でなくてはならないことがわか る．従って，劣化画像から劣化モデルを推定す るという操作が別個に必要となる． 同定問題は何らかの予備的知識の†に推定を 行なう．ここではぼけの種類は一様水平ぼけ， 雑音の種類は白色ガウス雑音と仮定し，ぼけの 大きさを推定することにより，ぼけのインパル ス応答を推定する． 2．2．2 一様水平ぼけ カメラなどで被写体を撮影する際に，得た画 像は手ぶれによって流れてしまうことがある． これが手ぶれによる劣化画像である．一様水平 ぼけと仮定しているので，ぼけの大きさ（次数） を推定することにより，インパルス応答が決定 可能である．例えば，ぼけの次数が2次の場合， 合，H・＝↓ 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0M次の場合，HM＝諸了
なる． 0 0 1 0 0 0 0 1 0 ，ぼけの次数が 0 ．． 0 0 0 03 一様水平ぼけの大きさ推定
3．1 一様水平ぼけの周波数特性
と 一様水平ぼけの周波数特性は以下のようにし て求められる． 2次の一様水平ぼけの場合，インパルス応答はH2−
（1） である．式（1）の2次元フーリエ変換を求める と以下のようになる． H（〆・画一；・一ゴW2（1＋・−」・・＋・一・」・・）（2） 4次の一様水平ぼけの場合は， 1H4＝
T
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 を2次元フーリエ変換して （3） H（e」 Wi，eJω2）−i・一・ゴ・・（1＋・一…＋・一・… ＋・−3」ω・＋e−4ゴω・）（4） となる．式（2）（4）より，M次のぼけHMの周波 数特性Hは式（5）のようになる． 昨・・，己㎏）＝M≒1・一￥ゴ姥（1＋・一ゴ軌 十．．．十e−（M−1）ごω1十e−Mゴω1） （5） 式（5）の，ω1方向の成分は，初項1，公比e一ゴω・の 等比数列である．よって，式（5）は以下のように書き直せる． H（ei7ω1，e」ω2） ニ？li lfit e− ￥」w2（1＿e−（M＋1）」ω・ 1＿．e−」ω1） 一砲ピ・・誓ω・一獅誓ω・） ・（1−c・s（M＋1）ω1＋元sin（M＋1）ω1 1−cosω1 −1−jsinω1）（6） よってM次の一様水平ぼけHの電力スペクト ルIH（ω1，ω2）12は，次式になる． IH（ω、，ω，）12＝ （M＋1）2 1 1−c・s（M＋1）ω1 （7） 1−cosω1 ω1≠0，2πで，かつ1−cos（M十1）ω1＝0と なるω1の時，IH（ω1，ω2）12はゼロになる．よっ て，111（ω1，ω2）12のゼロ点の数と，ぼけの大きさ （次数）は図一3，4，5の様に比例する．
3．2 周波数領域の劣化モデル2）
劣化モデルの行列演算式g＝Hf＋nは，周 波数領域では，次式で表される． G（ω1，ω2）ニH（ω1，ω2）F（ω1，ω2） ＋」V（ω1，ω2） （8） ここで，G（ω1，ω2），F（ω1，ω2），N（ω1，ω2）は，そ れぞれ，劣化画像，原画像，雑音の周波数特性 である．よって電力スペクトルによる劣化モデ ルは次式で表すことができる． IG（ω、，ω，）12＝IH（ω、，ω，）121F（ω、，ω，）12 −←IN（ω1，ω2）i2 （9） 雑音は平均値0の白色ガウス雑音と仮定してい るので，雑音の電力スペクトルIN（ω1，ω2）12の レベルは一定である．3．3 次数推定の原理
曇＿ 9．… trequency 図一3 2次の一様水平ぼけの周波数特性 Fig．3 Frequency Characteristics of uniform horizontal blur of order 2 曇＿ ≧ 図一4 4次の一様水平ぼけの周波数特性 Fig．4 Frequency Characteristics of unif（）rm horizontal blur of order 4 曇＿ 竃… freguetscツ 図一5 6次の一様水平ぼけの周波数特性 Fig．5 Frequency Characteristics of unifbrm horizontal blur of order 6 雑音が加わっていない場合，劣化画像の電力 スペクトルはぼけと原画像の電力スペクトルの 積となる． IG（ω、，ω2）12＝IH（ω、，ω，）121F（ω、，ω、）12（10） IH（ω1，ω2）12に次数に比例してゼロ点が存在す るので，IG（ω1，ω2）12にも同じ周波数で同じ数だ けゼロ点が存在する．雑音の加わった画像の場 合は，雑音を平均値0の白色ガウス雑音と仮定 しているので，ゼロ点が雑音の電力スペクトル IN（ω1，ω2）12分だけ持ち上がっている．よって， 劣化画像の電力スペクトルから雑音の電力スペ クトルを引くことによりゼロ点が検出できる． ぼけの電力スペクトルのゼロ点は次数に比例し て増加するので，劣化画像の電力スペクトルの ゼロ点の数を数えることによりぼけの次数が推 定できる．3．4 一様水平ぼけの次数推定の問題点
実際の劣化画像の電力スペクトルは，劣化画 像を2次元FFTして求めた離散系列であるた め，離散周波数点と一致しないゼロ点は明確な ゼロ点としては現れない．しかし，その場合で もゼロ点近傍の離散周波数にはゼロ付近まで落 ち込む点が存在する．図一・6は4次のぼけの劣化§ き 昆 iOOOO ioe O．i O．Oi frequency（ω1） 図一6 雑音の加わってない劣化画像の電力ス ペクトル（mandrill） Fig．6 Power spectrum of blured mandrill lmage 画像の電力スペクトルである．0一πの離散周波 数で2箇所落ち込む場所を確認できる．以降， ゼロ点とゼロ付近まで落ち込む点をまとめて落 下点と呼ぶことにする．以降，劣化画像の電力 スペクトルの落下点を検出することでぼけ次数 を推定する．その際，いくつかの間題が現れる． その問題と解決策を以下に示す． のレベルが異なりうまく閾値が決めれなくな る．そこでIL（ω1，ω2）12＝1 一 cos（ω1）の高域 通過特性のラプラシアンフィルタを劣化画像の 電力スペクトルに掛けることによって，レベル をほぼ一定とする．図一7よりラプラシアンフィ ルタを掛けた電力スペクトルはレベルが揃って いるのがわかる． 3．4．2 最大値フィルタ 電力スペクトルには細かな振動が無数にあり 落下点の検出が困難になる．そこで平滑化フィ ルタを電力スペクトルに掛け，細かな振動を取 り除き滑らかにする．しかし，線形の平滑化フィ ルタでは細かな振動は取れるが，急激な変化が 保存されず，本来の落下点までも消えてしまう． そこでエッジ保存性に優れる非線形の3近傍の 最大値フィルタ4）を繰り返し掛けることで細か な変動を取り除き，落下点の検出を容易にする． 図一8は電力スペクトルに3近傍の最大値フィル タを3回掛けたものである．細かな振動が取り 除け，急激な変化が保存されているのがわかる． 3．4．1 ラプラシアンフィルタ 原画像の電力スペクトルは，一般に低域に集 中しており，ぼけの周波数特性は低域通過特性 のため，劣化画像の電力スペクトルは図一6のよ うに低域が高域に比べて非常に高い．よって落 下点の検出の際に不都合が生じる．例えば，ス ペクトルを平滑化する時，最も低域に存在する 落下点は低域のレベルの高さが原因で消えてし まうことがある．また，閾値を決める時，落下点 §・・ § 三・ ftequency（（D i） 図一7 ラプラシアンフィルタを掛けた電力ス ペクトル Fig．7 Power spectrum mtering by Laplacian §・・ 皇1 f官equency（（D‘） 図一8最大値フィルタを掛けた電力スペクトル Fig．8 Power spectrum mtering by maximum 3．4．3 閾値の設定 図一8の電力スペクトルを見ると，最大値フィ ルタを掛けても，スペクトルの山の部分の細か なくぼみは消えにくいことがわかる．よって間 違った落下点を検出してしまい，次数を間違っ て推定してしまうおそれがある．そこで閾値を 決めて，山の部分の不要なくぼみをカットする． 図一9のように，閾値と電力スペクトルの交点の 数が最大となるレベルの一番低い値を閾値とす
る．ラプラシアンフィルタにより，望みの落下 点のレベルはほぼ一定になっているので，落下 点は保存されスペクトルの山の部分の細かなく ぼみを取り除くことができる．よって推定精度 が向上する． § 茎 § 』 閾値 freq cy
不要なくぼみ 望みの落下点
図一9 閾値 Fig．9 thrsuhold 3．4．4 落下点検出 劣化画像の電力スペクトルにはゼロ点が存在 しないものと仮定し，ぼけの周波数特性の落下 点と劣化画像の電力スペクトルの落下点の位置 と数が一致すれば，ぼけの次数を推定したこと とする． 原画像の電力スペクトルは高域のレベルが低 く，ぼけの周波数特性は低域通過特性である．し たがって劣化画像の電力スペクトルの高域のレ ベルは低く落下点を検出しにくい．また雑音が 乗った場合には高域の落下点は雑音の電力スペ クトルに埋もれてしまい検出できなくなる．そ こで，比較する落下点は低中域の落下点とし，．M 次の場合には，低域から数えて（誓＋1）／2個 の落下点を比較することとする．劣化画像の電 力スペクトルの落下点は最大値フィルタを掛け た回数分だけ幅を持たせ，高次のぼけの周波数 特性の落下点から順に比較していく．落下点の 数と周波数が一致すれば，周波数特性の次数を 劣化画像のぼけの次数とする．3．5 次数推定のアルゴリズム
前節の解決策をすべて組み込んだ実際の一様 水平ぼけの次数推定アルゴリズムは以下のよう になる． 1．劣化画像を2次元FFTする． 2．ω1方向の電力スペクトルを取り出す． 3．雑音の電力スペクトルを引きラプラシアン フィルタを掛ける． 4．3近傍の最大値フィルタを掛ける． 5．閾値を求めて，閾値以上の値を一定とする． 6．落下点を比較して，一致すれば推定次数と し，一致しなければ4−6を繰り返す． 以上より一様水平ぼけの次数は推定される．4 シミュレーション
4．1 シミュレーションの条件
使用する画像は256×256画素，8bit／pixel静 止画の白黒画像mandrill，lenna（図一10，12）であ る．ぼけの大きさが，2，4，6，8，10，12次の一様水平 ぼけを用いた． 観測画像には，信号対雑音比（SNR）が30dB になるように，白色ガウス雑音を加えた．ただ し，画像のSNRは次式で表される． variance of blurred image ［dB］SNR＝1010910
variance of the noise 雑音の電力スペクトルを求める際，雑音の分 散値を既知とした場合の他に，Wai Ho Punら の提案した雑音推定法3）を用いて雑音分散値を 推定した場合についても結果を示す．4．2 シミュレーション結果
表一1，2にmandrillのぼけの次数の推定結 果を，表一3，4にlennaのぼけの次数の推定結 果をそれぞれ示す．次に，図一10，12にman− drill，lennaの原画像を，図一11，13に，ぼけの次 数が8次，SNR＝ 30dBの観測画像をそれぞれ 示す．4．3 シミュレーション結果の考察
前節のシミュレーション結果より，ぼけの次 数はほぼうまく推定できていることがわかる． また，雑音の分散は推定結果にあまり影響がな いことがわかる． また，ω1方向の電力スペクトルを取り出す 時，ω2＝▲πの時が最も正確に推定されてい ることがわかる．これは，2次元の劣化画像の電 力スペクトルの高域部分のレベルは低域に比べ て低く，落下点が雑音の電力スペクトルに埋も れているためである．特に，lennaのような滑ら表一1mandrillの次数推定結果（雑音分散既知） Table−1 Result of estimated impulse re． sponse size of degraded mandrill image（Vari− ance of noise is already known） z の｝｝ 2 4 6 8 10 12 日の 1．31 1．22 1．59 1．12 1．10 1．05 ω2：＝π 堰ﾟil藪i 2222 4444 6666 8888 10 P0 W10 12 P2 P2 P2 表一2mandrillの次数推定結果（雑音分散推定） Table−2 Result of estimated impulse re− sponse size of degraded mandrill image（Vari− ance of noise is estimated） の↓1 2 4 6 8 10 12 同の 8．87 6．85 5．86 4．83 4．27 3．74 ω2＝π 奄撃堰ﾟ藪i 2222 4444 6666 8888 10 P0 P0 P0 12 P2 P2
m
N
よ 、つここ 刀く 表一3 1ennaの次数推定結果（雑音分散既知） Table−3 Result of estimated impulse re− sponse size of degraded lenna image（Variance of noise is already known） ノ の｝｝ 2 4 6 8 10 12 の 2．29 2．20 2．11 2．04 1．97 L91 ω2＝7「遠ｵi
8422 10 S44 6666 2888 841010 6212NN
よ、っ」こ 亦 表一4 1ennaの次数推定結果（雑音分散推定） Table−4 Result of estimated impulse re− sponse size of degraded lenna image（Variance of noise is estimated） ノ の｝1 2 4 6 8 10 12 の 2．06 1．56 1．31 1．18 1．09 LO1 ω2＝π蛾
8422 4444 10 U66 4888 16 Q1010 82124 N｝ で ヱ 、つここ 刀s かでコントラストのはっきりした画像の場合， 高域の電力スペクトルが低く，同じSNRでも 雑音の分散はmandrillなどと比べて大きくな るのでμ2＝πのとき，落下点が雑音に埋もれ てしまい，うまく推定できないことがわかる． また，ω2＝計の時は電力スペクトルの低域の レベルが高すぎて，ラプラシアンフィルタの効 果が低くなり，うまく推定できなくなる．よっ て，ω2＝れ付近の電力スペクトルを用いて推 定を行なうのが適当なことがわかる． また，表一1のように，ω2＝れで，雑音の分散 を真値とした時でも10次のぼけを間違って8 次と推定している．ここでは，雑音を平均値0 の白色ガウス雑音と仮定しているので，電力ス ペクトルは一定であるとみなしている．ところ が実際に計算に用いることができるのは電力ス ペクトルでなく，1つの標本系列のフーリエ変換 の2乗値であり，これは真値を境にして上下に 激しく振動している．よって，劣化画像の電力 スペクトルから雑音の電力スペクトルを引いた 場合に，落下点がなくなり，検出できなくなった と考えられる．5 おわりに
本研究では，劣化画像のスペクトルのゼロ点 を用いて，ぼけの大きさを推定し，インパルス 応答を決定する方法を提案した．ぼけを一様水 平ぼけと仮定することによって，ぼけの大きさ が推定できることをシミュレーションによって 示した．本方法では，水平方向の一様ぼけしか 推定できないため，斜めぼけなどにも対応した 推定方法を実現することが今後の課題である．参考文献
1）ReginaldL．Lagendijk，A．Murat Tekalp and Jan Biemond：“Maxdmum like− lihood image and blur identi丘cation：aunifying approach，，，OPTICAL
ENGINEERING，Vo1．29，No．5，
pp．422−435，May 1990． 2）Naoki Ono and Ryuzo Takiyama： “Estimation of Noise in Uniform M（卜 tion Blurred lmages for Restoration，， ，Proceedings of IEEE InternationalSymposium on Infbrmation Theory
and Its Applications，Vo1．1，pp．108− 112，1996． 3）Wai Ho Pun and Brian D．Jeffs： “Adaptive Image Restoration Using Model for Unknown Noise，，， IEEE Transactions on Image Processing， pp．1451−1456，1995． 4）棟安実治，雛元孝夫：“非線形デジタ ルフィルタ”，計測と制御，Vo1．35，No．9， pp．704−710，1996．図一10 原画像（mandrill） Fig．10 0riginal mandrill image． 図一12 原画像（1enna） Fig．12 0riginal lenna image． 図一11 劣化画像（mandrill 8次の一様水平ぼ けSNR＝30dB） Fig．11 Mandrill image degraded by a 1×9 unif（）rm motion blur with Gaussian Noise （SNR＝30dB）． 図一13 劣化画像（lenna 8次の一様水平ぼけ SNR＝30dB） Fig．13 Lenna image degraded by a 1×9 unif（）rm motion blur with Gaussian Noise （SNR＝30dB）．
