講義（山本裕樹）

この文書は読取専用です。文書のソースを閲覧することは可能ですが、変更はできません。もし変更したい場合は管理者に連絡してください。
====== 【数値情報処理b】第3回の答え ======
===== 課題 =====

gaussianfilter.R を書き換えて、$\sigma=1$ でサイズが $5\times 5$ のガウシアン・フィルタの畳み込み演算を5回くり返すプログラムを作成しなさい。

  * $5\times 5$ のフィルタを作成して ''correlate()'' を使うこと。
  * フィルタ作成には二項分布 $\mathrm{Bi}(4,0.5)$ を利用すること。

----
==== 解答例 ====

$\mathrm{Bi}(4,0.5)$ の確率密度関数は以下になる。
\[
f(x)=_{4}C_{x}(0.5)^{4}=_{4}C_{x}\times\frac{1}{16}
\]
$x=0,1,2,3,4$ の確率密度関数をベクトルの要素にすると
\[
\left(
\begin{array}{c}
_{4}C_{0}(0.5)^{4}\\
_{4}C_{1}(0.5)^{4}\\
_{4}C_{2}(0.5)^{4}\\
_{4}C_{3}(0.5)^{4}\\
_{4}C_{4}(0.5)^{4}
\end{array}
\right)
=
\left(
\begin{array}{c}
1/16\\
4/16\\
6/16\\
4/16\\
1/16
\end{array}
\right)
\]
このベクトルと、このベクトルを転置したものをかけあわせてフィルターを作る。
\begin{eqnarray*}
&&
\left(
\begin{array}{ccccc}
f(1,1)&f(2,1)&f(3,1)&f(4,1)&f(5,1)\\
f(1,2)&f(2,2)&f(3,2)&f(4,2)&f(5,2)\\
f(1,3)&f(2,3)&f(3,3)&f(4,3)&f(5,3)\\
f(1,4)&f(2,4)&f(3,4)&f(4,4)&f(5,4)\\
f(1,5)&f(2,5)&f(3,5)&f(4,5)&f(5,5)
\end{array}
\right)
\\
&=&
\left(
\begin{array}{c}
1/16\\
4/16\\
6/16\\
4/16\\
1/16
\end{array}
\right)
\left(
\begin{array}{ccccc}
1/16&4/16&6/16&4/16&1/16
\end{array}
\right)
\\
&=&
\frac{1}{256}
\left(
\begin{array}{ccccc}
1&4&6&4&1\\
4&16&24&16&4\\
6&24&36&24&6\\
4&16&24&16&4\\
1&4&6&4&1
\end{array}
\right)
\end{eqnarray*}

これが $\sigma=\sqrt{4\times 0.5\times(1-0.5)}=1$ の正規分布を近似したフィルタとして使える。
全体を 256 で割らなくても Imager では同じ効果になる。

<file rsplus report3.R [enable_line_numbers="true",highlight_lines_extra="12,13,14,15,16"]>
# 5x5のガウシアン・フィルタ

library(imager)

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

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

# ガウシアン・フィルタの畳み込み演算を繰り返す
f.gauss <- as.cimg(c(1,  4,  6,  4, 1,
                     4, 16, 24, 16, 4,
                     6, 24, 36, 24, 6,
                     4, 16, 24, 16, 4,
                     1,  4,  6,  4, 1), x = 5, y = 5)
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 ではフィルタを 256 で割っても結果は変わらない。