====== 【数値情報処理b】第3回 画像の平滑化 ======

===== 前回の課題 =====

[[ans2]]
===== ノイズと平滑化 =====

==== ノイズ ====
夜や暗い場所においてデジカメやスマホで写真を撮ると、画像にざらざらした邪魔なものがたくさん見えたことはないだろうか。
この邪魔なものは__**ノイズ（雑音）**__という。

下の写真はインドネシアのバリ島でデジカメで撮影した南十字星（サザンクロス）の写真である。
星が写っていない暗いところは真っ暗ではなく、拡大してよく見るとざらざらしたものが見える。

{{ bali_cross.png?nolink |バリ島で撮影した南十字星の写真}}

このようなざらざらと見えるものは、ほとんどが場所や大きさが不規則で予測できないノイズで、__**ランダムノイズ**__という。
ランダムノイズの原因としては、光の揺らぎによるものやデジカメ内部の電子のランダムな運動によるものがある。

----

==== 平滑化とコンポジット ====
ランダムノイズは場所的にランダムである。
ノイズの画素値とその周りのノイズではないところの画素値では急激な差がある。
この性質を利用してノイズを画像から取り除く手法を__**平滑化**__といい、
ボックス・フィルタ、メディアン・フィルタ、ガウシアン・フィルタなどの方法がある。
平滑化では、画像をぼかしてノイズを目立たなくする効果がある。

また、ランダムノイズは時間的にもランダムである。
動かない同じ被写体を何枚か撮影した場合、それぞれの画像の同じ場所にどのようなノイズが現れるかはランダムである。
この性質を利用して複数枚の画像を加算平均してノイズを小さくすることを__**コンポジット**__という。

平滑化は1枚の画像でできるが、コンポジットは同じ被写体を複数枚撮影しなければならないのがネックである。
しかし、コンポジットでは画像の枚数を増やすほどノイズを小さくできるというメリットがある。


===== GIMP による平滑化とコンポジット =====

GIMP で平滑化を行うには、ぼかしを使う。
メニューの ''**[フィルター]**''%%→%%''**[ぼかし]**'' で以下の中から選択する。

  * ガウスぼかし
  * 平均曲率ぼかし
  * メディアンぼかし

また、GIMP でコンポジットを行うには、複数枚の画像をレイヤーとして置き、レイヤーを「スクリーン」のモードで重ねることでできる。

===== ノイズ生成 =====
平滑化の練習に使用するために、ランダムノイズが入った画像を作ってみよう。

imager では ''imnoise()'' という関数で人工的にホワイトノイズと呼ばれるランダムノイズを生成できる。
<code rsplus>
> imnoise(x, y)
</code>
x は横幅、y は縦幅である。

次のプログラム whitenoise.R はグレースケール化した boats にホワイトノイズを加えている。
<file rsplus whitenosie.R[enable_line_numbers="true"]>
# ノイズ生成

library(imager)

# グレースケール画像
im <- grayscale(boats)

# ホワイトノイズを加えたノイズ画像
im.noise1 <- im + 0.1 * imnoise(width(im), height(im))
im.noise2 <- im + 0.1 * imnoise(width(im), height(im))

# 同時に3枚並べるレイアウト
layout(t(1 : 3))

plot(im, interpolate = FALSE, main = "元の画像")
plot(im.noise1, interpolate = FALSE, main = "ノイズ画像1")
plot(im.noise2, interpolate = FALSE, main = "ノイズ画像2")

# レイアウトを元に戻す
layout(1)
</file>

同じコマンドでノイズ画像を2枚生成しており、2枚を見比べるとノイズがランダムであることが分かる。

下の画像は whitenoise.R で生成した画像の同じ場所の部分を切り取った画像（上段）とその画像の $y=10$ の行の画素値のグラフ（下段）である。
元の画像の画素値が平坦なところは、ノイズ画像ではギザギザになっていることが分かる。

{{ whitenoise.png?nolink |ノイズ画像の比較}}

===== ボックス・フィルタ（移動平均法） =====
ある画素の周辺の画素の画素値の平均でその画素の画素値を置き換える手法を__**ボックス・フィルタ（移動平均法）**__という。
この方法は画像をぼかすことで細かいノイズもぼかして目立たなくしてしまう効果がある。
その反面、画像が全体的にぼやけて輪郭がはっきりしなくなってしまう。

例えば、下図の左側のような<fc #ff0000>ノイズを含む画素（中心の画素）</fc>を囲む全部で 9 つの画素を考える。
ノイズを含む画素は周辺の 8 つの画素とは大きく異なる画素値を持っている。

{{ boxfilter.png?nolink |ボックス・フィルタでノイズを消去}}

9 つの画素値の平均を計算すると
\[
\frac{1}{9}\times(0.1+0.2+0.1+0.4+0.9+0.3+0.1+0.3+0.3)=0.3
\]
である。
そこで、中心の画素値 0.9 を平均の 0.3 で置き換える。
置き換えた画素値は周辺の画素値の平均なので、置き換えても違和感なくノイズを消すことができる。

座標 $(x,y)$ の画素の周辺 $3\times 3$ の範囲の画素値の平均は以下である。

\begin{align*}
&\frac{1}{9}\{\ g(x-1,y-1)&&+g(x,y-1)&&+g(x+1,y-1)\\
&\ \ \ +g(x-1,y\ \ \ \ \ \ \ )&&+g(x,y\ \ \ \ \ \ \ )&&+g(x+1,y\ \ \ \ \ \ \ )\\
&\ \ \ +g(x-1,y+1)&&+g(x,y+1)&&+g(x+1,y+1)\}
\end{align*}

これは[[lecture2#畳み込み演算（コンボリューション）]]で表すことができる。
フィルタは以下のようになる。

\[
\left(
\begin{array}{ccc}
f(1,1)&f(2,1)&f(3,1)\\
f(1,2)&f(2,2)&f(3,2)\\
f(1,3)&f(2,3)&f(3,3)
\end{array}
\right)
=
\left(
\begin{array}{ccc}
1/9&1/9&1/9\\
1/9&1/9&1/9\\
1/9&1/9&1/9
\end{array}
\right)
=\frac{1}{9}
\left(
\begin{array}{ccc}
1&1&1\\
1&1&1\\
1&1&1
\end{array}
\right)
\]

次のプログラム boxfilter.R は boats をグレースケールのノイズ画像にし、$3\times 3$ のボックス・フィルタの畳み込み演算を5回繰り返している（1回では効果が小さいため）。

<file rsplus boxfilter.R[enable_line_numbers="true"]>
# ボックス・フィルタ

library(imager)

# グレースケール画像
im <- grayscale(boats)

# ノイズ画像
im.noise <- im + 0.1 * imnoise(width(im), height(im))

# ボックス・フィルタの畳み込み演算を5回繰り返す
f.box <- as.cimg(rep(1, 9), x = 3, y = 3) / 9
im.box <- im.noise
for (i in 1 : 5) {
  im.box <- correlate(im.box, f.box)
}

# 同時に3枚並べるレイアウト
layout(t(1 : 3))

plot(im, interpolate = FALSE, main = "元の画像")
plot(im.noise, interpolate = FALSE, main = "ノイズ画像")
plot(im.box, interpolate = FALSE, main = "ボックス・フィルタ")

# レイアウトを元に戻す
layout(1)
</file>

===== 小テスト① =====

[[KMS>|Moodle Server（非公式）]]で第3回の小テスト①を受験しなさい。

===== メディアン・フィルタ =====
ある画素の周辺の画素の画素値の中央値（メディアン）でその画素の画素値を置き換える手法を__**メディアン・フィルタ**__という。
ボックス・フィルタと比べて**<fc #ff0000>輪郭が損なわれにくい</fc>**という特徴がある。

例えば、下図の左側のようなノイズ（中心の画素）を囲む全部で 9 つの画素を考える。
ノイズは周辺の 8 つの画素とは大きく異なる画素値を持っている。

{{ medianfilter1.png?nolink |メディアン・フィルタでノイズを消去}}

9つの画素値は小さい順に並べると以下のようになる。

^         ^       ^       ^       ^       ^  中央値  ^       ^       ^       ^  ノイズ                ^
^ 順位    ^  1    ^  2    ^  3    ^  4    ^  5       ^  6    ^  7    ^  8    ^  9                     ^
^ 画素値  |  0.1  |  0.1  |  0.1  |  0.2  |  0.3     |  0.3  |  0.3  |  0.4  |  <fc #ff0000>0.9</fc>  |

中央値は 0.3 である。
ノイズは一番大きい（もしくは一番小さい）ので中央値にはなりにくい。
あとは、中心の画素値 0.9 を中央値 0.3 で置き換えるとノイズは消える。

輪郭については、例えば下図のような輪郭（青色の値）がある場合を考える。
輪郭はノイズと違い、周りの画素と連続的な画素値を持つ。

{{ medianfilter2.png?nolink |メディアン・フィルタで輪郭を処理}}

9 つの画素値は小さい順に並べると以下のようになる。

^         ^       ^       ^       ^                        ^  中央値                ^                        ^                        ^                        ^                        ^
^ 順位    ^  1    ^  2    ^  3    ^  4                     ^  5                     ^  6                     ^  7                     ^  8                     ^  9                     ^
^ 画素値  |  0.1  |  0.1  |  0.2  |  <fc #0000ff>0.8</fc>  |  <fc #0000ff>0.8</fc>  |  <fc #0000ff>0.8</fc>  |  <fc #0000ff>0.8</fc>  |  <fc #0000ff>0.9</fc>  |  <fc #0000ff>0.9</fc>  |

中央値は 0.8 となり、中心の画素値 0.9 を中央値 0.8 で置き換えてもあまり変わらない。

一方、平均は 0.6 なので、ボックス・フィルタでは輪郭に少し影響が出る。

メディアン･フィルタは畳み込み演算ではできないので、後ほど関数を用いたプログラムで行う。

===== 小テスト② =====

[[KMS>|Moodle Server（非公式）]]で第3回の小テスト②を受験しなさい。

===== ガウシアン・フィルタ =====
ある画素の周辺の画素の画素値に正規分布（ガウス分布）で重みをかけて足し合わせ、その画素の画素値を置き換える手法を__**ガウシアン・フィルタ**__という。
ボックス・フィルタと似たような効果になるが、ボックス・フィルタより自然な感じになる。

2変量正規分布（平均と共分散が $0$ で $\sigma_{x}=\sigma_{y}=\sigma$ の場合）の確率密度関数は以下である。
\[
f(x,y)=\frac{1}{2\pi\sigma^{2}}\exp\left(-\frac{x^{2}+y^{2}}{2\sigma^{2}}\right)
=\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{x^{2}}{2\sigma^{2}}\right)
\times\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{y^{2}}{2\sigma^{2}}\right)
\]
$\sigma$ の値が大きいほど分布は広がるので平滑化の効果は大きくなる。

{{ gauss2d.png?nolink |2変量正規分布}}

ガウシアン・フィルタは実数なので平滑化の計算に時間がかかる。
そこで計算を速くするために二項分布で整数に近似したフィルタの畳み込み演算がよく使われる。

__**ド・モアブル=ラプラスの定理**__によれば、
平均 $np$、標準偏差 $\sqrt{np(1-p)}$ の二項分布 $\mathrm{Bi}(n,p)$ は $n$ が大きいほど正規分布 $\mathrm{N}(np,np(1-p))$ に近づくという性質がある。
$\mathrm{Bi}(n,p)$ の確率密度関数 $f(x)$ は以下である。
\[
f(x)=_{n}C_{x}\;p^{x}(1-p)^{n-x}
\]
ここで $_{n}C_{x}$ は二項係数である。
$x$ は $x=0,1,2,\cdots,n$ の整数値をとる。


例えば、$3\times 3$ のフィルタを作る場合は $\mathrm{Bi}(2,0.5)$ を利用するとよい。
$\mathrm{Bi}(2,0.5)$ の確率密度関数は以下になる。
\[
f(x)=_{2}C_{x}(0.5)^{2}=_{2}C_{x}\times \frac{1}{4}
\]
$x=0,1,2$ の確率密度関数をベクトルの要素にすると
\[
\left(
\begin{array}{c}
_{2}C_{0}(0.5)^{2}\\
_{2}C_{1}(0.5)^{2}\\
_{2}C_{2}(0.5)^{2}
\end{array}
\right)
=
\left(
\begin{array}{c}
1/4\\2/4\\1/4
\end{array}
\right)
\]
このベクトルと、このベクトルを転置したものをかけあわせてフィルターを作る。
\begin{eqnarray*}
\left(
\begin{array}{ccc}
f(1,1)&f(2,1)&f(3,1)\\
f(1,2)&f(2,2)&f(3,2)\\
f(1,3)&f(2,3)&f(3,3)
\end{array}
\right)
&=&
\left(
\begin{array}{c}
1/4\\2/4\\1/4
\end{array}
\right)
\left(
\begin{array}{ccc}
1/4&2/4&1/4
\end{array}
\right)
\\
&=&
\left(
\begin{array}{ccc}
1/16&2/16&1/16\\2/16&4/16&2/16\\1/16&2/16&1/16
\end{array}
\right)
\\
&=&
\frac{1}{16}
\left(
\begin{array}{ccc}
1&2&1\\2&4&2\\1&2&1
\end{array}
\right)
\end{eqnarray*}
$1/16$ を除いた部分は整数のフィルタになっている。
これが $\sigma=\sqrt{2\times 0.5\times(1-0.5)}=\sqrt{0.5}\simeq 0.71$ の正規分布を近似したフィルタとして使える。

次のプログラム gaussianfilter.R は boats をグレースケールのノイズ画像にし、$3\times 3$ のガウシアン・フィルタの畳み込み演算を5回繰り返している。

<file rsplus gaussianfilter.R[enable_line_numbers="true"]>
# ガウシアン・フィルタ

library(imager)

# グレースケール画像
im <- grayscale(boats)

# ノイズ画像
im.noise <- im + 0.1 * imnoise(width(im), height(im))

# ガウシアン・フィルタの畳み込み演算を5回繰り返す
# Bi(2,0.5) を利用
f.gauss <- as.cimg(c(1, 2, 1, 2, 4, 2, 1, 2, 1), x = 3, y = 3) / 16
im.gauss <- im.noise
for (i in 1 : 5) {
  im.gauss <- correlate(im.gauss, f.gauss)
}

# 同時に3枚並べるレイアウト
layout(t(1 : 3))

plot(im, interpolate = FALSE, main = "元の画像")
plot(im.noise, interpolate = FALSE, main = "ノイズ画像")
plot(im.gauss, interpolate = FALSE, main = "ガウシアン・フィルタ")

# レイアウトを元に戻す
layout(1)
</file>



===== 平滑化の関数 =====
imager には、あらかじめ平滑化用の関数が用意されている。
パラメータをすぐに変えられるので便利である。

^  関数                           ^  機能                                                                         ^
| ''boxblur(cimg変数, boxsize)''  | ボックス・フィルタによる平滑化を行う。''boxsize'' はボックスのサイズである。  |
| ''medianblur(cimg変数, n)''     | メディアン・フィルタによる平滑化を行う。''n'' はフィルタのサイズである。      |
| ''isoblur(cimg変数, sigma)''    | ガウシアン・フィルタによる平滑化を行う。''sigma'' は標準偏差である。          |

----
==== 例 ====

smoothing.R は boats をノイズ画像にして平滑化を5回くり返すプログラムである。

<file rsplus smoothing.R[enable_line_numbers="true"]>
# 平滑化

library(imager)

# グレースケール画像
im <- grayscale(boats)

# ノイズ画像
im.noise <- im + 0.1 * imnoise(width(im), height(im))

# 平滑化（ボックス、メディアン、ガウシアン）を5回繰り返す
im.box <- im.noise
im.med <- im.noise
im.gauss <- im.noise
for (i in 1 : 5) {
  im.box <- boxblur(im.box, 3)
  im.med <- medianblur(im.med, 3)
  im.gauss <- isoblur(im.gauss, sqrt(0.5))
}

# 同時に3枚並べるレイアウト
layout(t(1 : 3))

plot(im.box, interpolate = FALSE, main = "ボックス・フィルタ")
plot(im.med, interpolate = FALSE, main = "メディアン・フィルタ")
plot(im.gauss, interpolate = FALSE, main = "ガウシアン・フィルタ")

# レイアウトを元に戻す
layout(1)
</file>


===== コンポジット =====
コンポジットは動かない同じ被写体を複数回撮影し、それらの画像を加算平均してノイズを小さくする。
<fc #ff0000>**被写体が動かないことが前提**</fc>となるので、天体の画像のノイズ除去に用いられることが多い。

次のプログラム composite.R は boats のノイズ画像を 9 回生成して加算平均するプログラムである。
関数 ''add()'' を使ってノイズ画像を1枚1枚加算して最後に平均している。
複数枚の画像を一度に平均する関数 ''average()'' もあるが、メモリがかなり必要になる。

<file rsplus composite.R[enable_line_numbers="true"]>
# コンポジット

library(imager)

# グレースケール画像
im <- grayscale(boats)

# 0で埋めたcimg変数
im.c <- imfill(dim = dim(im))

# 9枚のノイズ画像を加算
n <- 9
for (i in 1 : n) {
  # ノイズ画像
  im.noise <- im + 0.1 * imnoise(width(im), height(im))

  iml <- imlist(im.c, im.noise)
  im.c <- add(iml)
}
im.c <- im.c / n

# 同時に3枚並べるレイアウト
layout(t(1 : 3))

plot(im, interpolate = FALSE, main = "元の画像")
plot(im.noise, interpolate = FALSE, main = "最後のノイズ画像")
plot(im.c, interpolate = FALSE, main = "コンポジット")

# レイアウトを元に戻す
layout(1)
</file>

{{ composite.png?nolink |9枚の画像のコンポジットを画素値で比較 }}

元の画像と比べると 9 枚でもノイズをかなり除去できており、元の画像をある程度再現できることが分かる。

===== 課題 =====
[[KMS>|Moodle Server（非公式）]] で第3回の課題を行いなさい。

締め切り：2025年7月2日（水）20時