2004 年度卒業論文　手ぶれ劣化画像の復元における復元フィルタ評価

(1)

2004 年度卒業論文

手ぶれ劣化画像の復元における復元フィルタ評価

2005 年 2 月 2 日

指導教授：大石進一教授

早稲田大学理工学部情報学科学籍番号： 1G01P038-0

栗原崇

(2)

第 1 _{章序論}

1.1 背景

近年、画像工学の分野において、劣化した画像から原画像を復元する研究が多くなされている。これは防犯や医療などの分野での応用が期待されている。撮影された画像には、様々な原因によって歪みやぼけが生じる。例えば、ビデオカメラなどで撮影された画像は、カメラレンズの焦点のずれ、手ぶれ、シャッタースピード、撮影素子の性能などの影響から、実物と比較して画質が劣化することは避けられない。画像復元とはこのように劣化した画像からできるだけ原画像に近い画像を復元することである。

劣化した画像は一般的に原画像と点分布関数

(PSF)

の畳み込み

(コンボ

リューション)のモデルによって得られる。PSFが既知であれば、逆フィルタや最小二乗フィルタ

(ウィナーフィルタ)

などのフィルタリング処理方法を用いて原画像を復元できることが知られている。このような問題は、逆畳み込み

(デコンボリューション)

と呼ばれる。

1.2 本論文の目的

1.1

節で述べた通り、PSFが既知であれば、復元フィルタリングすることで原画像を復元できる。しかし、

PSF

にもそれぞれ劣化特性があり、それに対する有効な復元フィルタリングが存在するはずである。そこで本論文は、劣化画像復元方法についてまとめ、復元画像の収束性及び画質の改善を行うことを目的とする。

本論文では

2

種類の制限範囲を設けた逆フィルタと最小二乗フィルタ

(ウィナーフィルタ)

を

C

言語で実装し、手ぶれ画像に限定した各々の劣

化画像復元方法の有効性について検証する。また画像の輝度値量子化の際に生じる誤差を考慮した復元方法についても検証していく。

(6)

1.3 本論文の構成

本論文の構成は以下の通りである。

第

2

章では、画像復元理論として、劣化画像、逆フィルタ、最小二乗

フィルタ

(ウィナーフィルタ)

について述べる。

第

3

章では、手ぶれ画像復元の準備として、実際に

PSF

を求め、対応する逆フィルタ、最小二乗フィルタ

(ウィナーフィルタ)

を構成する。

第

4

章では、画像の輝度値量子化の際に生じた誤差を考慮した復元方法について述べる。

第

5

章では、画像復元の実行結果を示し、検証を行う。

第

6

章では、結論について述べる。

尚、付録として本論文で作成したプログラムリストを示す。

5

(7)

第 2 _{章画像復元の理論}

2.1 はじめに

本章では、画像復元の理論となる劣化画像モデル、逆フィルタ、最小二乗フィルタ

(ウィナーフィルタ)

について述べる。

2.2 劣化画像のモデル

イメージングシステム

(対象物を観測する装置や画像を生成するシステ

ム)によって得られる画像

g(x, y )

は、次の式で表される。

g(x, y) =

Z _∞

−∞

Z _∞

−∞

h(x, y, α, β, γ)f (α, β)dαdβ (2.1)

ただし、f(α, β)は対象物を表す原画像、h(·,

·, ·, ·)

は撮影条件なども含めたイメージングシステムの特性を表す点広がり関数

(PSF)

である。イメージングシステムは線形で位置不変な特性を持つものとすると、式

(2.1)

は

2

次元の線形たたみこみ演算

g(x, y) =

Z _∞

−∞

Z _∞

−∞

h(x − α, y − β)f (α, β)dαdβ (2.2)

となり、これを

g(x, y) = f(x, y) ∗ h(x, y) (2.3)

と表現する。

また、式

(2.3)

の両辺に

2

次元フーリエ変換を施すことによって、空間

周波数領域では、

G(µ, ν) = F{g(x, y)}

= F{f (x, y) ∗ h(x, y)}

= H(µ, ν )F (µ, ν) (2.4)

(8)

と表すことができる。ここで、µ, νはそれぞれ

x, y

方向への空間周波数を表している。また、H(µ, ν)を光学的伝達関数

(OTF)

という。

カメラなどで観測する際には雑音が含まれるはずであるので、その観

測雑音を

n(x, y)、そのフーリエ変換を N (µ, ν)

とすると、劣化の過程を

示す式

(2.2)

及び式

(2.4)

は、

g(x, y) =

Z _∞

−∞

Z _∞

−∞

h(x − α, y − β)f(α, β)dαdβ + n(x, y) (2.5) G(µ, ν) = H(µ, ν)F (µ, ν ) + N (µ, ν) (2.6)

となる。

ディジタル処理を行う場合には、式

(2.2)

はディジタル畳み込み演算となり、原画像及び観測画像は辞書式配列によって並べ替えられてそれぞれベクトル

f

及び

g

となり、次のように書き換えられる。

g = [H]f + n (2.7)

ただし、nは観測雑音やディジタル化を行うときに生じた誤差等を含めた雑音ベクトルである。

2.3 復元画像のモデル

復元画像

ˆf

は、観測画像

g

に復元フィルタ

[B ]

を作用させてその出力として得られる。復元フィルタをアナログ領域、または周波数領域で考えると次のようになる。

復元フィルタの点広がり関数を

b(x, y)、そのフーリエ変換を B(µ, ν)

とすると、復元画像

f ˆ (x, y)

とそのフーリエ変換

F ˆ (µ, ν)

は、

f ˆ (x, y) =

Z _∞

−∞

Z _∞

−∞

b(x − α, y − β)g(α, β)dαdβ (2.8) F ˆ (µ, ν ) = [H(µ, ν)F (µ, ν ) + N (µ, ν)]B(µ, ν) (2.9)

で得られる。式

(2.8)

の処理を式

(2.5)

に対して逆畳み込み

(deconvolution)

という。

2.3.1 逆フィルタ

復元フィルタの伝達関数

B (µ, ν)

を

B(µ, ν) = 1

H(µ, ν ) (2.10)

7

(9)

とすると、復元画像

F ˆ (µ, ν )

は、

F ˆ (µ, ν ) = F (µ, ν) + N (µ, ν)

H(µ, ν) (2.11)

となる。式

(2.10)

の特性をもつ復元フィルタを

H(µ, ν )

の逆フィルタ

(in-

verse filter)

という。雑音が存在しなければ、逆フィルタによって原画像

F (µ, ν )

は完全に復元される。

しかし雑音の存在は避けられず、H(µ, ν)が零かまたは零に非常に近い値をとる周波数領域があると、わずかの雑音でも式

(2.11)

の第

2

項が著しく拡大されて復元画像は大きく乱れてしまう。

ディジタル処理においても、式

(2.7)

で逆フィルタに相当する

[H]

の逆行列

[H]

⁻¹が存在すれば、復元画像

ˆf

は次式で得られる。

ˆf = f + [H]

⁻¹

(2.12)

しかし、一般に

[H]

⁻¹は存在しないことが多く、仮に存在しても正確に求めることは極めて困難になる。このような状況でなんらかの復元処理を強行すると式

(2.12)

の右辺第

2

項が拡大されて、著しく劣化した復元画像しか得られなくなる。

このような状況に陥る問題を悪条件問題

(ill-conditioned problem)

または不良設定問題

(ill-posed problem)

といい、画像復元はほとんど全てが悪条件問題となる。悪条件問題で、雑音の拡大を抑える方策を正則化

(reguralization)

という。従って、画像復元ではいかにして正則化処理を

繰り込むかが大きなポイントとなる。

2.3.2 最小二乗フィルタ ( ウィナーフィルタ )

復元フィルタで得られた画像における復元誤差を測る評価基準として、

次の式で定義される原画像と復元画像との間の平均

2

乗誤差を用いる場合を考える。

Ek²k

²

= E[{f(x, y) − f ˆ (x, y )}

²

] (2.13)

ただし、E[·]は平均操作を示す。観測画像に作用させてこの平均

2

乗誤差を最小にする画像を得る復元フィルタを最小

2

乗フィルタ

(least square

filter)

またはウィナーフィルタ

(Wiener filter)

という。

(10)

原画像や付随する雑音成分が弱定常場に属すると仮定できるときには、

最小

2

乗フィルタの空間周波数特性は次式で与えられる。

B (µ, ν ) = H

^∗

(µ, ν)

|H(µ, ν)|

²

+ S

_n

(µ, ν)/S

_f

(µ, ν) (2.14)

ただし、∗は複素共役を表し、S_f

(µ, ν)

及び

S

_n

(µ, ν)

は原画像及び雑音のパワースペクトル密度で、画像及び雑音の自己相関関数

R

f

(x, y)

及び

R

_n

(x, y)

のフーリエ変換として次のように定められる。

S

_f

(µ, ν ) =

Z _∞

−∞

Z _∞

−∞

R

_f

(x, y) exp{−j2π(µx + νy)}dxdy (2.15) S

_n

(µ, ν ) =

Z _∞

−∞

Z _∞

−∞

R

_n

(x, y) exp{−j 2π(µx + νy)}dxdy (2.16)

また、これらが未知の場合、式

(2.14)

の近似として、

B(µ, ν) ' H

^∗

(µ, ν)

|H(µ, ν)|

²

+ Γ (2.17)

を用いる。ここで、Γは定数である。

式

(2.7)

で得られるディジタル画像の場合、最小

2

乗フィルタは次のよ

うにして求められる。復元フィルタ

[B ]

によって復元画像

ˆf

は、次式で求められる。

ˆf = [B]g (2.18)

原画像との差

²

の

2

乗平均値、すなわち平均

2

乗誤差を最小にする

[B]

を求める。すなわち、

min E[²²

^T

] = min E[(f − ˆf)(f − ˆf)

^T

] (2.19)

となる

[B]

を求める。

式

(2.7),(2.18)

を用いて、

E[(f − [B]g)(f − [B]g)

^T

]

= E[{ff

^T

− [B ]([H]ff

^T

+ nf

^T

) − (ff

^T

[H] + fn

^T

)[B]

+[B]([H]ff

^T

[H]

^T

+ nf [H]

^T

+ [H]fn

^T

+ nn

^T

)[B]

^T

}] (2.20)

加わる雑音は画像と無相関であると考えてよいので、

E[fn

^T

] = E[n

^T

f ] = 0 (2.21)

9

(11)

となり、f 及び

n

の相関行列を

E [ff

^T

] = [R

_{f f}

], E[nn

^T

] = [R

_nn

] (2.22)

と表すと、式

(2.20)

は次のようになる。

E[²²

^T

] = [R

_{f f}

] − 2[B][H][R

_{f f}

] + [B][H][R

_{f f}

][H][B]

^T

+ [B ][R

_nn

][B]

^T

(2.23)

そこで、

∂E [²²

^T

]

∂ [B] = 0 (2.24)

から、

[B] = [R

_{f f}

][H]

^T

([H][R

_{f f}

][H]

^T

+ [R

_nn

])

⁻¹

(2.25)

が求める最小

2

乗フィルタとなる。原画像ベクトル

f

は、2次元画像を辞書式配列に従って並べかえて作られたものであるから、fの相関行列は次のようになる。

[R

_{f f}

] =

¯¯

R

₀₀

R

₀₁

. . . R

_0(M−1)

R

₁₀

R

₁₁

. . . ...

... ... ... ...

R

_(M₋₁₎₀

. . . R

_(M−1)(M₋₁₎

¯¯

(2.26)

式

(2.26)

を構成する各行列

R

ij は

N × N

行列で、次式で表される。

R

_ij

= E[f

_i

f

_j^T

] (2.27)

ただし、fiと

f

jは

i

番目と

j

番目の行ベクトルである。Rij

= R

jiとなる

ので式

(2.26)

は対称行列となる。

雑音が存在しない場合には式

(2.14)

及び式

(2.25)

のフィルタは、いずれも逆フィルタと一致する。最小

2

乗フィルタでは、式

(2.13)

や式

(2.19)

の評価基準から出発することによって復元フィルタの中に雑音の成分が繰り込まれて自動的に正則化が行われている。例えば、式

(2.14)

におい

て

H(µ, ν) = 0

となるような領域でも、B(µ, ν)となって、逆フィルタで

生じたように雑音が極端に増幅されて演算が不安定になり復元画像が大きく乱れることはない。

(12)

2.4 むすび

本章では、画像復元の理論となる劣化画像モデル、逆フィルタ、最小二乗フィルタ

(ウィナーフィルタ)

について述べた。

次章では、手ぶれ画像復元の準備として、実際に

PSF

(ウィナーフィルタ)

を構成する。

11

(13)

第 3 _{章手ぶれ画像の復元}

3.1 はじめに

本研究では、現実の劣化画像モデルとして、手ぶれ画像に焦点をあてる。そこで本章では、手ぶれ画像復元の準備として、実際に

PSF

を求め、

対応する逆フィルタ、最小二乗フィルタ

(ウィナーフィルタ)

を構成する。

3.2 手ぶれ劣化画像の生成

写真を撮る際に生じるカメラの手ぶれによるぼけを考える。対象とな

る風景

f(x, y)

とカメラとの相対運動が、x方向、y方向にそれぞれ

α(t)、

β(t)

の移動で表されるとし、シャッターの開いている時間を

T

とすれば、

フィルム面上に得られる画像

g(x, y)

は、

g(x, y) =

Z ^T

2

−^T₂

f(x − α(t), y − β(t))dt (3.1)

と表すことができる。この式の両辺をフーリエ変換すると、

G(µ, ν)

=

Z Z ÃZ ^T

2

−^T₂

f(x − α(t), y − β(t))dt

!

exp{−j2π(µx + νy)}dxdy

=

Z ^T

2

−^T₂

dt

Z Z

f (x − α(t), y − β(t)) exp{−j 2π(µx + νy)}dxdy (3.2)

ここで、

u = x − α(t), v = y − β(t) (3.3)

と置くと、

G(µ, ν ) =

Z ^T

2

−^T₂

dt

Z Z

f (u, v) exp{−j 2π(µ(u + α(t)) + ν(v + β(t))}dudv

(14)

=

Z ^T

2

−^T₂

dt

Z Z

f (u, v) exp{−j 2π(µu + νv)} exp{−j 2π(µα(t) + νβ(t))}dudv

=

Z ^T

2

−^T₂

exp{−j 2π(α(t)µ + β(t)ν)}dt × F (µ, ν) (3.4)

ここで、

H(µ, ν) =

Z ^T

2

−^T₂

exp{−j2π(α(t)µ + β(t)ν)}dt (3.5)

と置くと、

G(µ, ν) = H(µ, ν )F (µ, ν) (3.6)

となるので、式

(3.5)

は手ぶれによるぼけの伝達関数だということがわかる。

本研究では、手ぶれは等速直線運動とする。V を動きの速度、

θ

を動きの方向とすると、

α(t) = V cos θ × t, β(t) = V sin θ × t (3.7)

と置けて、

H(µ, ν) =

Z ^T

2

−^T₂

exp{−j2π(V cos θ × tµ + V sin θ × tν}dt

= sin πV (µ cos θ + ν sin θ)T πV (µ cos θ + ν sin θ)

= T sin c{πV (µ cos θ + ν sin θ)T } (3.8)

以上、手ぶれによるぼけの伝達関数を求めることができた。

3.3 逆フィルタによる手ぶれ劣化画像の復元

式

(2.10)、式 (3.8)

より求める復元フィルタの伝達関数

B(µ, ν )

は、

B (µ, ν) = 1 H(µ, ν)

= 1

T sin c{πV (µ cos θ + ν sin θ)T } (3.9)

となる。この復元フィルタを用いることで、復元画像

F ˆ (µ, ν)

が以下のように求まる。

F ˆ (µ, ν ) = G(µ, ν)B(µ, ν)

13

(15)

= G(µ, ν)

T sin c{πV (µ cos θ + ν sin θ)T } (3.10)

これを逆フーリエ変換することで復元画像

f ˆ

が得られる。

しかし、式

(3.8)

より

H(µ, ν) = 0

となる点があるので式

(2.11)

より画像を正確に復元できない。そこで、制限を加えることでこの問題を回避する。

3.3.1 領域制限を加えた逆フィルタリング

逆フィルタを適用する領域を制限してやることで回避する。ある

r

に対して、

B(µ, ν) =





1

H(µ,ν)

µ

²

+ ν

²

≤ r

²

1 µ

²

+ ν

²

> r

²

(3.11)

3.3.2 無効制限範囲を加えた逆フィルタリング

逆フィルタを無効にする

H(µ, ν)

の制限範囲を設けることで回避する。

ある

²

に対して、

B (µ, ν) =





1

H(µ,ν)

|H(µ, ν)| > ²

0 |H(µ, ν)| ≤ ² (3.12)

3.4 最小二乗フィルタ ( ウィナーフィルタ ) による手ぶれ劣化画像の復元

式

(2.14)、式 (3.8)

より求める復元フィルタの伝達関数

B(µ, ν )

は、

B(µ, ν) = H

^∗

(µ, ν )

|H(µ, ν)|

²

+ S

n

(µ, ν)/S

f

(µ, ν)

= T sin c{πV (µ cos θ + ν sin θ)T }

|T sin c{πV (µ cos θ + ν sin θ)T }|

²

+ S

n

(µ, ν )/S

f

(µ, ν) ' T sin c{πV (µ cos θ + ν sin θ)T }

|T sin c{πV (µ cos θ + ν sin θ)T }|

²

+ Γ (3.13)

(16)

となる。この復元フィルタを用いることで、復元画像

F ˆ (µ, ν)

が以下のように求まる。

F ˆ (µ, ν ) = G(µ, ν)B(µ, ν)

= G(µ, ν)T sin c{πV (µ cos θ + ν sin θ)T }

|T sin c{πV (µ cos θ + ν sin θ)T }|

²

+ S

_n

(µ, ν)/S

_f

(µ, ν) ' G(µ, ν)T sin c{πV (µ cos θ + ν sin θ)T }

|T sin c{πV (µ cos θ + ν sin θ)T }|

²

+ Γ (3.14)

これを逆フーリエ変換することで復元画像

f ˆ

が得られる。

3.5 むすび

本章では、手ぶれ劣化に対応する

PSF

を求め、それを復元する逆フィルタ、最小二乗フィルタ

(ウィナーフィルタ)

を構成した。

次章では、量子化誤差を縮小する復元方法について述べる。

15

(17)

第 4 章量子化誤差縮小復元

4.1 はじめに

ディジタル画像を作る際、標本化によって空間的に離散化された画素に分解するが、各画素の値

(濃度値,

輝度値)に関してはアナログのままである。そこで、連続した濃淡値を離散的な整数値に変換する操作が量子化である。そこには当然誤差が伴うわけで、この誤差が画像復元におい

ては式

(2.5)、式 (2.11)

のように影響してくるわけである。本章では、こ

の量子化誤差を縮小する復元方法について述べる。

4.2 _{量子化誤差縮小復元}

量子化前の劣化画像

g(x, y) ˜

は、

˜

g(x, y) = F

⁻¹

{F{f (x, y)}H(µ, ν)} (4.1)

量子化後の劣化画像

g(x, y)

は、

g(x, y) = 5˜ g(x, y)

= ˜ g(x, y) + n(x, y) (4.2)

ここで、5は下向きの丸め操作、n(x, y)は量子化誤差を表す。

以上のことから、以下のことが導き出せる。

0 ≤ g(x, y) ≤ g(x, y) ˜ ≤ g(x, y) + 1 ≤ (最大階調値)

⇐⇒ F

⁻¹

{F{g(x, y)}B(µ, ν)} ≤ F

⁻¹

{F{˜ g(x, y)}B(µ, ν)}

≤ F

⁻¹

{F{g(x, y) + 1}B (µ, ν )}

or F

⁻¹

{F{g(x, y)}B(µ, ν)} ≥ F

⁻¹

{F{˜ g(x, y)}B(µ, ν)}

≥ F

⁻¹

{F{g(x, y) + 1}B (µ, ν )}

⇐⇒ F

⁻¹

{F{g(x, y)}B(µ, ν)} ≤ f (x, y )

(18)

≤ F

⁻¹

{F{g(x, y) + 1}B (µ, ν )}

or F

⁻¹

{F{g(x, y)}B(µ, ν)} ≥ f (x, y )

≥ F

⁻¹

{F{g(x, y) + 1}B (µ, ν )} (4.3)

このことから、

f ˆ (x, y) = 5 F

⁻¹

{F{g(x, y)}B(µ, ν)} + F

⁻¹

{F{g(x, y) + 1}B(µ, ν)}

2 (4.4)

とすることで、量子化誤差を縮小した復元画像

f(x, y) ˆ

が得られる。

4.3 _むすび

本章では、量子化誤差を縮小する復元方法について述べた。

次章では、手ぶれ劣化画像の復元結果を示す。

17

(19)

第 5 _{章実行結果}

5.1 はじめに

本章では、実際に手ぶれ劣化画像を作成し、各々の復元方法を用いて復元した結果を示す。また、見た目では捉えられない精度を測る尺度として

PSNR(Peak Signal to Noise Ratio)

を導入し、定量的に精度比較を行う。

5.2 _実行環境

CPU : Celeron 2.00GHz

Memory : 768MB

OS : Windows XP

(20)

5.3 画像フォーマット

本研究で扱う画像ファイルは、白黒階調画像に対するシンプルかつ汎用性の高く、理解し易い画像フォーマットである、pgm(portable graymap)

形式

(バイナリ)

のファイルを採用する。

5.4 原画像

本研究で利用する原画像を以下のような

256 × 256

の白黒

256

階調画像とする。

図

5.1:

原画像

19

(21)

5.5 劣化画像生成

式

(3.8)

のそれぞれ

T, V, θ

に値を与えることで手ぶれ伝達関数が求め

られ、式

(4.1)

を逆フーリエ変換することで劣化画像は得られる。以下は

PSF

とそれに対応する劣化画像の例である。尚、本研究では雑音を量子化誤差だけとする。

図

5.2: T = 0.1, V = 1, θ = 0

図

5.3: T = 0.1, V = 1, θ = 0

図

5.4: T = 0.1, V = 1, θ = π/2

図

5.5: T = 0.1, V = 1, θ = π/2

(22)

5.5.1 手ぶれの速度、角度を固定 (V = 1, θ = 0)

図

5.6: T = 0.01, V = 1, θ = 0

図

5.7: T = 0.01, V = 1, θ = 0

図

5.8: T = 0.1, V = 1, θ = 0

図

5.9: T = 0.1, V = 1, θ = 0

21

(23)

図

5.10: T = 1, V = 1, θ = 0

図

5.11: T = 1, V = 1, θ = 0

5.5.2 カメラの露光時間、手ぶれの角度を固定

(T = 0.1, θ = 0)

図

5.12: T = 0.1, V = 0.1, θ = 0

図

5.13: T = 0.1, V = 0.1, θ = 0

(24)

図

5.14: T = 0.1, V = 1, θ = 0

図

5.15: T = 0.1, V = 1, θ = 0

図

5.16: T = 0.1, V = 10, θ = 0

図

5.17: T = 0.1, V = 10, θ = 0

23

(25)

5.6 復元過程

手ぶれ伝達関数として、式

(3.8)

に

T = 0.1, V = 1, θ = 0

を適用した場合の復元画像を表示していく。

5.6.1 領域制限を加えた逆フィルタリング

図より、領域制限を加えた逆フィルタリングでは手ぶれ画像を復元することはできない。

図

5.18:

復元画像

(r = 5)

図

5.19:

復元

PSF(r = 5)

図

5.20:

復元画像

(r = 9)

図

5.21:

復元

PSF(r = 9)

(26)

図

5.22:

復元画像

(r = 10)

図

5.23:

復元

PSF(r = 10)

図

5.24:

復元画像

(制限なし)

図

5.25:

復元

PSF(制限なし)

25

(27)

5.6.2 無効制限範囲を加えた逆フィルタリング

図より、

² = 0.004

の時、最も良好な復元結果が得られることがわかる。

図

5.26:

復元画像

(² = 0.010)

図

5.27:

復元

PSF(² = 0.010)

図

5.28:

復元画像

(² = 0.004)

図

5.29:

復元

PSF(² = 0.004)

(28)

図

5.30:

復元画像

(² = 0.001)

図

5.31:

復元

PSF(² = 0.001)

図

5.32:

復元画像

(制限なし)

図

5.33:

復元

PSF(制限なし)

27

(29)

5.6.3 最小二乗フィルタ ( ウィナーフィルタ )

図より、Γ = 0.00010の時、最も良好な結果が得られるが、無効制限範囲を加えた逆フィルタリングと比べると劣る。

図

5.34:

復元画像

(Γ = 0.00100)

図

5.35:

復元

PSF(Γ = 0.00100)

図

5.36:

復元画像

(Γ = 0.00050)

図

5.37:

復元

PSF(Γ = 0.00050)

(30)

図

5.38:

復元画像

(Γ = 0.00010)

図

5.39:

復元

PSF(Γ = 0.00010)

図

5.40:

復元画像

(Γ = 0.00005)

図

5.41:

復元

PSF(Γ = 0.00005)

29

(31)

図

5.42:

復元画像

(Γ = 0.00001)

図

5.43:

復元

PSF(Γ = 0.00001)

(32)

5.6.4 量子化誤差縮小復元

復元アルゴリズムとして最も良好な復元結果が得られた、無効制限範囲を加えた逆フィルタリングを用いて、量子化誤差縮小復元の結果を示す。

図

5.44:

復元画像

(² = 0.010)

図

5.45:

復元

PSF(² = 0.010)

図

5.46:

復元画像

(² = 0.005)

図

5.47:

復元

PSF(² = 0.005)

31

(33)

図

5.48:

復元画像

(² = 0.004)

図

5.49:

復元

PSF(² = 0.004)

図

5.50:

復元画像

(² = 0.001)

図

5.51:

復元

PSF(² = 0.001)

(34)

5.7 復元結果

手ぶれ伝達関数を変化させた時の復元結果を、最も良好な復元結果が得られた、量子化誤差縮小した無効制限範囲を加えた逆フィルタリングによる復元画像を例にして以下に示す。

5.7.1 手ぶれの速度、角度を固定

(V = 1, θ = 0)

図

5.52:

劣化画像

(T = 0.05)

図

5.53:

復元画像

(T = 0.05)

図

5.54:

劣化画像

(T = 0.25)

図

5.55:

復元画像

(T = 0.25)

33

(35)

5.7.2 カメラの露光時間、手ぶれの角度を固定 (T = 0.1, θ = 0)

図

5.56:

劣化画像

(V = 0.5)

図

5.57:

復元画像

(V = 0.5)

図

5.58:

劣化画像

(V = 2.5)

図

5.59:

復元画像

(V = 2.5)

(36)

5.8 PSNR 比較

5.8.1 PSNR(Peak Signal to Noise Ratio)

原画像と復元画像との差異を評価する尺度として、PSNR(Peak Signal

to Noise Ratio)

を用いる。

P SNR = 10 log

₁₀

µ

(最大階調値)

²

sqrt({f(x, y) − f ˆ (x, y)}

²

)

¶

(dB) (5.1)

5.8.2 復元過程における量子化誤差縮小の効果

無効制限範囲を変化させていく復元過程において、量子化誤差縮小の効果を

PSNR

比較で検証する。

inverse2

inverse2 2

量子化誤差縮小した無効制限範囲を加えた逆フィルタリング

² inverse2 inverse2 2 0.001 26.852376 26.853394 0.002 33.251105 33.311163 0.003 33.320486 33.380404 0.004 33.340275 33.398870 0.005 33.339408 33.369434 0.010 33.280122 33.280385

表

5.1: PSNR(dB)

35

(37)

5.8.3 手ぶれ伝達関数に対する量子化誤差縮小復元の効果

手ぶれ伝達関数の変数である、カメラの露光時間、手ぶれの速度を変化させ、各々の復元結果の量子化誤差縮小の効果を

PSNR

比較で検証する。

inverse2

inverse2 2

量子化誤差縮小した無効制限範囲を加えた逆フィルタリング

T inverse2 inverse2 2

0.01 ∗ ∗

0.05 33.414531 33.443846 0.10 33.340275 33.398870 0.15 33.338891 33.368337 0.20 33.263328 33.292726 0.25 33.144230 33.170711

1.00 ∗ ∗

表

5.2: PSNR(dB)(V = 1, θ = 0)

V inverse2 inverse2 2

0.1 ∗ ∗

0.5 33.416530 33.447274 1.0 33.340275 33.398870 1.5 33.342936 33.372222 2.0 33.298598 33.326756 2.5 33.200220 33.227025

10.0 ∗ ∗

表

5.3: PSNR(dB)(T = 0.1, θ = 0)

(38)

5.9 むすび

本章では、実際に手ぶれ劣化画像を作成し、各々の復元方法を用いて復元した結果を示し、PSNR比較を行った。

次章では、この結果を元に考察する。

37

(39)

第 6 _{章結論}

6.1 統括

本論文では

2

種類の制限付き逆フィルタと最小二乗フィルタ

(ウィナー

フィルタ)を実装し、手ぶれ画像における各々の劣化画像復元方法の有効性について検証した。また、画像の輝度値量子化の際に生じる誤差を縮小する復元方法を提案し、PSNRの比較によりその有効性を実証した。

第

2

章では、画像復元理論として、劣化画像、逆フィルタ、最小二乗

フィルタ

(ウィナーフィルタ)

について述べた。

第

3

章では、手ぶれ画像復元の準備として、実際に

PSF

(ウィナーフィルタ)

を構成した。また、零点問題を回避するために制限付き逆フィルタを

2

種類構成した。

第

4

章では、画像の輝度値量子化の際に生じた誤差を縮小させる復元方法について述べた。

第

5

章では、画像復元の実行結果を示し、検証を行った。

本論文の構成は以上の通りであった。

6.2 考察

まず、領域制限を用いた逆フィルタでは、手ぶれによるぼけ画像を復元することはできないことがわかる。これは手ぶれ伝達関数が急激に

0

に収束していくのではなく、sinカーブを描きながら緩やかに収束していくことが原因と考えられる。一方、無効制限範囲を用いた逆フィルタでは、

手ぶれ伝達関数の値によらず、零点問題を回避し、良好な復元結果を得ることができる。また、最小二乗フィルタ

(ウィナーフィルタ)

はぼけ画像を復元することはできるが、無効制限範囲を用いた逆フィルタに比べると、見た目で違いがすぐにわかる程、精度に差が生じる。よって量子化誤差に対しては、最小二乗フィルタは有効な復元方法ではないことがわかる。

(40)

量子化誤差縮小の効果は見た目ではわからないが、PSNR(dB)値の増加と定量的結果として現れているので、原画像との誤差が小さい、高精度な画像復元ができているといえる。

また、手ぶれの度合いが大きい、つまりカメラの露光時間、手ぶれの速度が大きい場合、今回用いた何れの復元方法を用いても復元することはできなかった。これは手ぶれによるぼけが原画像全体に及んでしまい、

劣化画像に復元に関する原画像の情報が十分に保持されていなかったためと考えられる。

6.3 結論

手ぶれによるぼけ画像の復元には、最小二乗フィルタ

(ウィナーフィル

タ)よりも量子化誤差縮小付きの無効制限範囲を用いた逆フィルタの方が有効であることが示せた。また、手ぶれのぼけが全体に及んでしまうような場合、原画像の情報が十分に保持されないため、逆フィルタや最小二乗フィルタ

(ウィナーフィルタ)

で復元することは不可能であることが示せた。

39

(41)

謝辞

本研究を進めるにあたり、終始丁寧な御指導及び御激励を賜り、その他多くの面でも色々と御面倒を見て下さり、御助言を与えて下さいました大石進一教授に深く感謝致します。

また、終始丁寧な御指導と御教示をして下さいました大石研究室助手、

丸山晃佐氏、荻田武史氏、に大いに感謝致します。

また、日常生活において色々とお世話になりました、大石研究室博士課程

1

年、尾崎克久氏、並びに大石研究室修士課程

2

年、大山博毅氏、森山敦史氏、並びに大石研究室修士課程

1

年、坂内太郎氏、島本誠氏、関本竜平氏、徳永克久氏、福地健氏、細田幸裕氏、横山大輔氏、に深く感謝致します。

また、意見の交換や協力などをして下さいました大石研究室

4

年、有滝貴広氏、岩田揚平氏、大成高顕氏、川崎文裕氏、佐藤友規氏、平野晃氏、

ト詣各氏、水島裕氏、吉岡毅氏、吉田昂央氏、に深く感謝致します。

最後に、研究だけでなく日常の生活の中でお世話になりました大石研究室の皆様に深く感謝致します。

(42)

参考文献

1. N.Mastronardi, P.Lemmerling, and S.Van Huffel

：

Fast structured to- tal least squares algorithm for solving the basic deconvolution prob- lem,SIAM Journal on Matrix Analysis and Applications,22(2):533- 553,(2000)

2.

斎藤恒雄：画像処理アルゴリズム、近代科学社、1998年

3.

長尾真：ディジタル画像処理、近代科学社、1978年

4.

大石進一：フーリエ解析、岩波書店、2002年

41

(43)

付録 A _{プログラムリスト}

A.1 画像の入出力

A.1.1 pgm.h

/**********************************************************/

/* pgm.h */

/* pgm.c 用ヘッダファイル */

/**********************************************************/

/* 定数宣言 */

#define MAX_IMAGESIZE 1024 /* 想定する縦・横の最大画素数 */

#define MAX_BRIGHTNESS 255 /* 想定する最大階調値 */

#define GRAYLEVEL 256 /* 想定する階調数(=最大階調値+1) */

#define MAX_FILENAME 256 /* 想定するファイル名の最大長 */

#define MAX_BUFFERSIZE 256 /* 利用するバッファ最大長 */

/* 大域変数の宣言 */

/* 画像用配列１，画像用配列２ */

unsigned char image1[MAX_IMAGESIZE][MAX_IMAGESIZE], image2[MAX_IMAGESIZE][MAX_IMAGESIZE];

int x_size1, y_size1, /* image1 の横画素数，縦画素数 */

x_size2, y_size2; /* image2 の横画素数，縦画素数 */

/* 関数のプロトタイプ宣言 */

void load_image_data( ); /* 画像読み込み用関数 */

void save_image_data( ); /* 画像書き込み用関数 */

A.1.2 pgm.c

/***************************************************************/

/* pgm.c */

/* pgm ファイルの入出力 */

/***************************************************************/

#include <stdio.h>

#include <stdlib.h>

(44)

#include "pgm.h"

/***************************************************************/

/* pgm 画像，横画素数，縦画素数のデータをファイルから読み込み，*/

/* image1[ ][ ]，x_size1，y_size1 にそれぞれ代入する． */

/***************************************************************/

void load_image_data( ) {

char file_name[MAX_FILENAME]; /* ファイル名用の文字配列 */

char buffer[MAX_BUFFERSIZE]; /* データ読み込み用作業変数 */

FILE *fp; /* ファイルポインタ */

int max_gray; /* 最大階調値 */

int x, y; /* ループ変数 */

/* 入力ファイルのオープン */

printf("---\n");

printf(" モノクロ階調画像入力ルーチン\n");

printf("---\n");

printf("ファイル形式は pgm, バイナリ形式とします．\n");

printf("入力ファイル名 (*.pgm) : ");

scanf("%s",file_name);

fp = fopen( file_name, "rb" );

if ( NULL == fp ){

printf("その名前のファイルは存在しません．\n");

exit(1);

}

/* ファイルタイプ(=P5)の確認 */

fgets( buffer, MAX_BUFFERSIZE, fp );

if ( buffer[0] != ’P’ || buffer[1] != ’5’ ){

printf("ファイルのフォーマットが P5 とは異なります．\n");

exit(1);

}

/* x_size1, y_size1 の代入（#から始まるコメントは読み飛ばす） */

x_size1 = 0;

y_size1 = 0;

while ( x_size1 == 0 || y_size1 == 0 ){

if ( buffer[0] != ’#’ ){

sscanf( buffer, "%d %d", &x_size1, &y_size1 );

} }

/* max_gray の代入（#から始まるコメントは読み飛ばす） */

max_gray = 0;

while ( max_gray == 0 ){

if ( buffer[0] != ’#’ ){

sscanf( buffer, "%d", &max_gray );

} }

43

(45)

/* パラメータの画面への表示 */

printf("横の画素数 = %d, 縦の画素数 = %d\n", x_size1, y_size1 );

printf("最大階調値 = %d\n",max_gray);

if ( x_size1 > MAX_IMAGESIZE || y_size1 > MAX_IMAGESIZE ){

printf("想定値 %d x %d を超えています．\n", MAX_IMAGESIZE, MAX_IMAGESIZE);

printf("もう少し小さな画像を使って下さい．\n");

exit(1);

}

if ( max_gray != MAX_BRIGHTNESS ){

printf("最大階調値が不適切です．\n");

exit(1);

}

/* 画像データを読み込んで画像用配列に代入する */

for ( y = 0; y < y_size1; y ++ ){

for ( x = 0; x < x_size1; x ++ ){

image1[y][x] = (unsigned char)fgetc( fp );

} }

printf("データは正しく読み込まれました．\n");

printf("---\n");

fclose(fp);

}

/*****************************************************************/

/* image2[ ][ ], x_size2, y_size2 のデータを，それぞれ */

/* pgm 画像，横画素数，縦画素数としてファイルに保存する． */

/*****************************************************************/

void save_image_data( ) {

char file_name[MAX_FILENAME]; /* ファイル名用の文字配列 */

FILE *fp; /* ファイルポインタ */

int x, y; /* ループ変数 */

/* 出力ファイルのオープン */

printf("---\n");

printf(" モノクロ階調画像（pgm形式）出力ルーチン\n");

printf("---\n");

printf("出力ファイル名 (*.pgm) : ");

scanf("%s",file_name);

fp = fopen(file_name, "wb");

/* ファイル識別子 "P5" を先頭に出力する */

fputs( "P5\n", fp );

/* # で始まるコメント行（省略可能） */

fputs( "# Created by Image Processing\n", fp );

/* 画像の横幅，縦幅の出力 */

fprintf( fp, "%d %d\n", x_size2, y_size2 );

/* 最大階調値の出力 */

fprintf( fp, "%d\n", MAX_BRIGHTNESS );

(46)

/* 画像データの出力 */

for ( y = 0; y < y_size2; y ++ ){

for ( x = 0; x < x_size2; x ++ ){

fputc( image2[y][x], fp );

} }

printf("データは正しく出力されました．\n");

printf("---\n");

fclose(fp);

}

A.2 フーリエ変換

A.2.1 FFT1.h

/**********************************/

/* FFT1.h */

/* FFT1.c用ヘッダファイル */

/**********************************/

#define FFT_MAX 1024

#define PI 3.141592653589 /* 円周率 */

int calc_power_of_two( int number );

void make_initial_data( double *re_part, double *im_part, int num_of_data, int power );

void FFT1( double *re_part, double *im_part, int num_of_data, int flag );

A.2.2 FFT1.c

/***************************************/

/* FFT1.c */

/* 1次元 FFT */

/***************************************/

#include "FFT1.h"

/***************************************/

/* number が２の何乗であるか調べて返す */

/* 例）8 --> 3, 16 --> 4, 32 --> 5,... */

/***************************************/

int calc_power_of_two( int number ) {

int power;

power = 0;

45

(47)

while ( number != 1 ){

power ++;

number = number >> 1;

/* 右方向に１ビットシフト（２で割るのと同じ） */

}

return power;

}

/*********************************************************************/

/* データを並び替えて時間の間引きによる FFT の初期データを作る */

/* double *re_part, */

/* *im_part; 元データ(実数部，虚数部)の先頭アドレス */

/* int num_of_data; データ総数 */

/* int power; num_of_data は２の power 乗 */

/*********************************************************************/

void make_initial_data( double *re_part, double *im_part, int num_of_data, int power ) {

int i, j; /* 制御変数 */

int ptr, offset; /* 元データの要素番号を決めるための変数 */

int new_ptr; /* 計算後のデータを代入する配列の要素番号 */

int dft; /* DFT点数 ( i のループで半分ずつになる） */

double re_buf[FFT_MAX],

im_buf[FFT_MAX]; /* 計算結果用配列 */

dft = num_of_data;

for ( i = 1; i < power; i ++ ){

/* i は回数を制御する変数（計算には用いていない） */

new_ptr = 0; offset = 0;

while( new_ptr < num_of_data ){

ptr = 0;

while( ptr < dft ){

re_buf[new_ptr] = *( re_part + offset + ptr );

im_buf[new_ptr] = *( im_part + offset + ptr );

new_ptr ++;

ptr = ptr + 2;

if ( ptr == dft ) ptr = 1;

}

offset = offset + dft;

}

/* 計算結果を元のデータ配列にコピーする */

for ( j = 0; j < num_of_data; j ++ ){

*( re_part + j ) = re_buf[j];

*( im_part + j ) = im_buf[j];

}

dft = dft / 2;

} }

/*********************************************************************/

(48)

/* データの FFT (flag = 1), 逆 FFT (flag = -1) を行う関数 */

/* double *re_part, */

/* *im_part; 元データ(実数部，虚数部)の先頭アドレス */

/* int num_of_data, flag; データ総数，FFT・逆FFT を決めるフラグ */

/*********************************************************************/

void FFT1( double *re_part, double *im_part, int num_of_data, int flag ) {

int i,j,k,power; /* 制御変数，べき乗 */

double unit_angle,step_angle; /* 角度用の変数 */

int dft,half,num_of_dft; /* DFT点数，その半数，DFT 実行回数 */

int num_out,num_in1,num_in2; /* DFTの入出力信号線の番号 */

double re_buf,im_buf,angle; /* 実数部，虚数部，角度用作業変数 */

/* 次は計算結果を一時的に保存するための作業用配列 */

static double re_part_new[FFT_MAX],im_part_new[FFT_MAX];

/* 逆FFT ( flg = -1 ) のとき各データを num_of_data で割り */

/* 複素共役をとる */

if ( flag == -1){

for ( i = 0; i < num_of_data; i ++ ){

*( re_part + i ) = *( re_part + i ) / num_of_data;

*( im_part + i ) = - *( im_part + i ) / num_of_data;

} }

/* num_of_data が２の何乗(power)であるか調べる */

power = calc_power_of_two( num_of_data );

/* 初期データのための元データの順番の入れ替え */

make_initial_data( re_part, im_part, num_of_data, power );

/* 2点DFT, 4点DFT, 8点DFT, ... の順に実行する */

unit_angle = 2.0 * PI / num_of_data;

dft = 2; /* = DFT の点数の初期値，２倍してゆく */

for ( i = 0; i < power; i ++ ){

/* DFT 点の DFT を行う */

/* i=0 -> 2点, i=1 -> 4点, i=2 -> 8点, .... */

num_of_dft = num_of_data / dft; /* DFT の実行回数 */

step_angle = unit_angle * num_of_dft;

half = dft / 2;

for ( j = 0; j < num_of_dft; j ++ ){

angle = 0.0;

for ( k = 0; k < dft; k ++ ){

/* No.num_in1 と No.num_in2 から No.num_out を出力 */

/* （係数が１の方の信号を No.num_in1 としている） */

num_out = j * dft + k;

if ( k < half ){

num_in1 = num_out;

num_in2 = num_in1 + half;

} else {

num_in2 = num_out;

num_in1 = num_out - half;

47

(49)

}

/* 実数部(re_)・虚数部(im_)に分けて計算 */

re_buf = *( re_part + num_in2 );

im_buf = *( im_part + num_in2 );

re_part_new[num_out] = *( re_part + num_in1 ) + re_buf * cos(angle) + im_buf * sin(angle);

im_part_new[num_out] = *( im_part + num_in1 ) + im_buf * cos(angle) - re_buf * sin(angle);

/* 角度を更新 */

angle = angle + step_angle;

} }

/* 計算後のデータを元の配列に戻す */

*( re_part + j ) = re_part_new[j];

*( im_part + j ) = im_part_new[j];

}

dft = dft * 2; /* DFT の点数を２倍に更新 */

}

/* 逆FFT ( flg = -1 ) のとき各データの複素共役をとる */

if ( flag == -1 ){

*( im_part + j ) = - *( im_part + j );

} } }

A.2.3 FFT2.h

/**************************************/

/* FFT2.h */

/* FFT2.c用ヘッダファイル */

/**************************************/

/* 大域変数 */

double data[FFT_MAX][FFT_MAX]; /* 元データの実数部 */

double jdata[FFT_MAX][FFT_MAX]; /* 元データの虚数部 */

int num_of_data; /* データの要素番号の最大値 */

/* 関数のプロトタイプ */

void FFT2( int flag );

A.2.4 FFT2.c

/****************************/

/* FFT2.c */

(50)

/* 2次元 FFT */

/****************************/

#include "FFT1.h"

#include "FFT2.h"

/**********************************************************/

/* 元データを data, jdata, num_of_data に代入してから使う */

/* flag = 1 : ２次元 FFT, flag = -1 : ２次元逆 FFT */

/**********************************************************/

void FFT2( int flag ) {

int i, j; /* ループ変数 */

static double re[FFT_MAX], im[FFT_MAX]; /* 作業変数 */

re[j] = data[i][j];

im[j] = jdata[i][j];

}

FFT1( re, im, num_of_data, flag );

data[i][j] = re[j];

jdata[i][j] = im[j];

} }

re[j] = data[j][i];

im[j] = jdata[j][i];

}

FFT1( re, im, num_of_data, flag );

data[j][i] = re[j];

jdata[j][i] = im[j];

} } }

A.3 画像データ ⇐⇒ フーリエ変換用データ

A.3.1 convertdata.h

/*********************************/

/* convertdata.h */

/* convertdata.c用ヘッダファイル */

49

(51)

/*********************************/

/* 関数のプロトタイプ宣言 */

void make_original_data( ); /* 画像データをFFT用データへ変換 */

void make_output_image( ); /* FFT用データを画像データへ変換 */

A.3.2 convertdata.c

/******************************************/

/* convertdata.c */

/* 画像データ <--> FFT用データの切り替え */

/******************************************/

#include <stdio.h>

#include <stdlib.h>

#include "pgm.h"

#include "FFT1.h"

#include "FFT2.h"

/**************************************************************/

/* image1[y][x] のデータを data[y][x], jdata[y][x] に代入する */

/**************************************************************/

void make_original_data( ) {

if ( x_size1 != y_size1 ){

printf("縦と横の画素数が異なっているので終了します．\n");

exit(-1);

} else {

printf("\n読み込んだ画像を元データに直します．\n");

num_of_data = x_size1;

data[i][j] = (double)image1[i][j];

jdata[i][j] = 0.0;

} } } }

/*****************************************************************/

/* データ data[y][x], jdata[y][x] を image2[y][x] に直す */

/* 計算値の最小値, 最大値を求めて [0,255] になるように正規化する */

/*****************************************************************/

void make_output_image( )

(52)

{

int x, y; /* 制御変数 */

double max, min; /* 計算値の最小値, 最大値 */

max = 0.0;

min = MAX_BRIGHTNESS;

printf("\n逆FFT変換後の画像を作成します．\n");

x_size2 = num_of_data; y_size2 = num_of_data;

for ( y = 0; y < y_size2; y ++ ){

for ( x = 0; x < x_size2; x ++ ){

if ( data[y][x] < 0 ) data[y][x] = 0;

if ( data[y][x] < min ) min = data[y][x];

if ( data[y][x] > max ) max = data[y][x];

} }

for ( y = 0; y < y_size2; y ++ ){

for ( x = 0; x < x_size2; x ++ ){

data[y][x] = MAX_BRIGHTNESS/(max-min)*(data[y][x]-min);

image2[y][x] = (unsigned char)data[y][x];

} } }

A.4 劣化画像生成

A.4.1 blur.c

/*********************************************************************/

/* blur.c */

/* カメラと被写体との相対的な動きによるボケ画像生成 */

/*********************************************************************/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include "pgm.h"

#include "FFT1.h"

#include "FFT2.h"

#include "convertdata.h"

/************************/

/* sinc関数 */

/************************/

double sinc(double x){

double n;

if ( x == 0 ) n = 1;

51

(53)

else n = sin(x)/x;

return n;

}

/****************************************/

/* 周波数領域でぼけフィルタリングを行う */

/****************************************/

void blur( ) {

double T; /* カメラの露光時間 */

double V; /* 動きの速度 */

double x; /* 動きの方向 */

double h[num_of_data][num_of_data]; /* 伝達関数 */

printf("カメラの露光時間 T :");

scanf("%lf",&T);

printf("動きの速度 V :");

scanf("%lf",&V);

printf("動きの方向 θ :");

scanf("%lf",&x);

/* 伝達関数を求める */

h[i][j] = T*sinc(PI*V*T*(j*cos(x)+i*sin(x)));

} }

printf("\nFFT後の係数に対するぼけフィルタリングを行います．\n");

data[i][j] *= h[i][j];

jdata[i][j] *= h[i][j];

} } }

/************************/

/* メイン関数 */

/************************/

main( ) {

load_image_data( ); /* 原画像の読み込み */

make_original_data( ); /* 元のデータを作成する */

FFT2( 1 ); /* ２次元 FFT の実行 */

blur( ); /* 周波数領域でのぼけフィルタリング */

FFT2( -1 ); /* ２次元逆FFT の実行 */

make_output_image( ); /* 変換データを画像に直す */

(54)

save_image_data( ); /* 出力画像を保存する */

return 0;

}

A.5 復元画像生成

A.5.1 inverse.c

/************************************************************/

/* inverse.c */

/* 周波数領域制限を用いた逆フィルタによるカメラと被写体との */

/* 相対的な動きによるボケ画像復元 */

/************************************************************/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include "pgm.h"

#include "FFT1.h"

#include "FFT2.h"

#include "convertdata.h"

/************************/

/* sinc関数 */

/************************/

double sinc(double x){

double n;

if ( x == 0 ) n = 1;

else n = sin(x)/x;

return n;

}

/***************************************************************/

/* 周波数領域で逆フィルタリングを行う */

/* ただし、領域制限を設けて制限外では逆フィルタの値を 1 とする */

/***************************************************************/

void inverse( ) {

double T; /* カメラの露光時間 */

double V; /* 動きの速度 */

double x; /* 動きの方向 */

double r; /* 周波数領域制限半径 */

double b[num_of_data][num_of_data]; /* 逆フィルタ */

2004 年度卒業論文 手ぶれ劣化画像の復元における復元フィルタ評価