【数値情報処理b】第2回の答え

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

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

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

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

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

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

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

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

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

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

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

imager では imnoise() という関数で人工的にホワイトノイズと呼ばれるランダムノイズを生成できる。

> imnoise(x, y)

x は横幅、y は縦幅である。

次のプログラム whitenoise.R はグレースケール化した boats にホワイトノイズを加えている。

whitenosie.R

# ノイズ生成
 
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)

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

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

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

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

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*}

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

\[ \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回では効果が小さいため）。

boxfilter.R

# ボックス・フィルタ
 
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)

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

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

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

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

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

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

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

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

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

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$ の値が大きいほど分布は広がるので平滑化の効果は大きくなる。

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

ド・モアブル=ラプラスの定理によれば、 平均 $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回繰り返している。

gaussianfilter.R

# ガウシアン・フィルタ
 
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)

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

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

smoothing.R

# 平滑化
 
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)

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

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

composite.R

# コンポジット
 
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)

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

Moodle Server（非公式） で第3回の課題を行いなさい。

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