====== 【数値情報処理b】第5回 2値化と物体検出 ======

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

[[ans4]]


===== しきい値処理 =====
画像の各画素について、画素値が一定値以上の場合にそれを抜き出す方法がある。
この基準となる値を__**しきい値**__という。
具体的には、画像の画素値がしきい値未満なら画素値 0（黒）、しきい値以上なら画素値 1（白）を割り当てる処理を行う。
このように画素値を 0 か 1 のどちらかに割り当てる処理を__**2値化**__という。

{{ boats_binary.png?nolink |2値化}}

しきい値が 0.3 の場合の画素の2値化
|  <fc #ff0000>0.4</fc>  |  <fc #4682b4>0.2</fc>  |  <fc #ff0000>0.3</fc>  ^  %% %%\\ %%→%%\\ 2値化  |  <fc #ff0000>1</fc>  |  <fc #4682b4>0</fc>  |  <fc #ff0000>1</fc>  |
|  <fc #ff0000>0.5</fc>  |  <fc #4682b4>0.1</fc>  |  <fc #ff0000>0.8</fc>  | :::                     |  <fc #ff0000>1</fc>  |  <fc #4682b4>0</fc>  |  <fc #ff0000>1</fc>  |
|  <fc #4682b4>0.1</fc>  |  <fc #ff0000>0.7</fc>  |  <fc #ff0000>0.6</fc>  | :::                     |  <fc #4682b4>0</fc>  |  <fc #ff0000>1</fc>  |  <fc #ff0000>1</fc>  |

次のプログラム threshold.R は、画像 nameplate.jpg をしきい値 ''thr'' を使って2値化する。

<file rsplus threshold.R[enable_line_numbers="true"]>
# 2値化

library(imager)

# グレースケール画像
im <- grayscale(load.image("https://www.kitp.org/numericalipb/nameplate.jpg"))

# しきい値
thr <- 0.3

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

# imの画素値がしきい値以上なら画素値1とする
im.bin[im >= thr] <- 1

layout(t(1 : 2))
plot(im, main = "元の画像")
plot(im.bin, main = "2値化")
layout(1)
</file>

しきい値 ''thr'' を 0 から 1 の範囲で色々変えると物体（ここでは文字）が見えたり見えなかったりすることが分かる。
物体を見分けやすくするための最適なしきい値はどうやって決めればよいだろうか？

===== GIMPによる2値化 =====

GIMP で2値化を行うには、メニューの ''**[色]**''%%→%%''**[しきい値]**'' を選択する。 

グレースケール画像の場合は ''**[明度]**'' で表示されるヒストグラムにおいて、▲ のしきい値の位置を手動で変更したり、自動調整ができる。


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

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


===== 2値化の方法 =====

ここでは2値化のための適切なしきい値を求める方法をいくつか挙げる。

----
==== モード法 ====
画素値のヒストグラムを描き、その谷をしきい値とする方法を__**モード法**__という。
ヒストグラムとは、データの値の範囲を階級で区切り、それぞれの階級中に含まれるデータの個数を度数とし、階級値ごとに度数を柱状グラフで表したものである。

物体と背景で濃度が大きく違う場合に谷が出やすい。

{{ detection_mode.png?nolink |モード法}}

R では ''hist()'' を使えばお手軽にヒストグラムを描くことができる。

cimg 変数の im のヒストグラムを描いてみよう。
<code rsplus>
> hg <- hist(im, breaks = seq(0, 1, 1 / 256))
</code>
''seq(0, 1, 1 / 256)'' は 0 から 1 までを 1 / 256 の等幅で 256 分割する境界値を生成する。
引数 ''breaks'' に与えることで 0 から 1 までの範囲の画素値を 256 の階級に分割できる。

''hist()'' の戻り値には以下の要素が含まれる。

^  ラベル  ^  要素         ^
| ''breaks''      | 階級の境界値  |
| ''counts''      | 度数          |
| ''density''     | 確率密度（相対度数を階級幅で割ったもの）      |
| ''mids''        | 階級値        |

ヒストグラムがギザギザして分かりにくい場合は、度数分布の前後 2 つの階級の度数で__**移動平均**__をとってヒストグラムの平滑化を行うと良い。
<code rsplus>
> x <- filter(hg$counts, rep(1 / 5, 5))
> hg$counts <- replace(x, is.na(x), 0)
> plot(hg)
</code>
''hg$counts'' は ''hist()'' で得られた度数の数列である。
これを ''filter()'' を使って移動平均を計算して ''x'' に代入している。

<hidden 移動平均の計算の詳細>
''rep(1 / 5, 5)'' は 1 / 5 を 5 回繰り返した数を生成する。
''filter()'' では 5 つの 1 / 5 を hg$counts のそれぞれの度数の前後と併せてかけて和をとっている。

例えば、''hg$counts'' が
<code rsplus>
5   4   6   1   9   3   8
</code>
の移動平均を考えてみよう。
前から 3 つめの度数 6 では、前後の 2 つずつの度数との移動平均は 5.0 になる。
\[
\frac{5}{5}+\frac{4}{5}+\frac{6}{5}+\frac{1}{5}+\frac{9}{5}=\frac{5+4+6+1+9}{5}=5.0
\]
同様に前から 4 つめの度数 1 の前後の 2 つずつの度数との移動平均は 4.6 になる。
\[
\frac{4}{5}+\frac{6}{5}+\frac{1}{5}+\frac{9}{5}+\frac{3}{5}=\frac{4+6+1+9+3}{5}=4.6
\]

一方、一番前の度数 5 の前後の 2 つずつの度数との移動平均を計算しようとすると、度数 5 の前には度数が無いため、移動平均の結果は ''NA''（欠損値）になってしまう。

以上より、''x <- filter(hg$counts, rep(1 / 5, 5))'' の結果は
<code rsplus>
NA  NA  5.0  4.6  5.4  NA  NA
</code>
になる。

''NA'' だとヒストグラムが描けないので、''replace()'' で ''NA'' を ''0'' に置換している。
''is.na()'' は ''NA'' かどうかを判断する関数である。
''replace(x, is.na(x), 0)'' の結果は
<code rsplus>
0.0  0.0  5.0  4.6  5.4  0.0  0.0
</code>
になる。

</hidden>
\\

移動平均をとった度数を見てみよう。
<code rsplus>
> hg$counts
Time Series:
Start = 1 
End = 256 
Frequency = 1 
  [1]    0.0    0.0 1027.2  971.6  935.6  864.6  805.4  733.2  652.0  587.0
 [11]  499.2  418.8  338.8  267.8  214.4  179.2  151.0  138.2  123.8  118.0
 [21]  110.8  105.6   97.0   96.6   89.0   88.8   86.4   85.4   87.0   88.0
 [31]   84.2   82.4   83.0   74.8   72.8   73.2   75.0   72.6   77.2   76.6
 [41]   77.2   75.6   76.0   72.2   73.6   69.0   68.4   68.8   68.8   66.8
 [51]   68.8   69.8   70.6   73.6   76.8   83.6   86.2   88.0   92.2   89.6
 [61]   86.8   88.6   93.0   93.4   97.6   99.8  103.0  100.2   98.8  103.2
 [71]  107.2  111.4  115.0  124.6  128.0  131.0  137.0  147.8  155.2  172.6
 [81]  205.0  243.0  302.2  405.6  595.8  976.4 1582.4 2305.4 2953.8 3382.2
 [91] 3455.8 3144.6 2617.8 2027.8 1534.4 1151.8  896.0  698.6  581.2  487.2
[101]  415.6  364.4  317.6  275.2  230.4  194.8  158.0  128.8  101.8   78.6
[111]   61.4   43.4   30.6   20.8   13.6    7.2    5.0    3.0    1.8    1.4
[121]    1.0    0.4    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0
[131]    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0
以下省略
</code>
谷の底は 50 番目付近なので、その階級値を見てみると
<code rsplus>
> hg$mids[50]
[1] 0.1933594
</code>
これより、しきい値は 0.19 あたりにすると良さそうである。

threshold.R でしきい値 ''thr'' を 0.19 にして実行してみよう。

----
==== p-タイル法 ====
物体が背景に比べて明るい場合、画像全体に対する背景の面積（画素数）の割合があらかじめ $p$ であると分かっているとする。
ヒストグラムで暗い方から $100p$ % の割合となるところをしきい値とする方法を__**p-タイル法**__という。
このときのしきい値はいわゆる __**$100p$ パーセンタイル**__である。

{{ detection_p-tile.png?nolink |p-タイル法}}


次のプログラム p-tile.R は画像を割合 p の p-タイル法で2値化する。

<file rsplus p-tile.R[enable_line_numbers="true",highlight_lines_extra="12"]>
# p-タイル法

library(imager)

# グレースケール画像
im <- grayscale(load.image("https://www.kitp.org/numericalipb/nameplate.jpg"))

# 割合
p <- 0.3

# しきい値の決定
thr <- quantile(im, p)
print(thr)

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

# imの画素値がしきい値以上なら画素値1とする
im.bin[im >= thr] <- 1

layout(t(1 : 2))
plot(im, main = "元の画像")
plot(im.bin, main = "2値化")
layout(1)
</file>

12 行目にある ''quantile(im, p)'' は引数 ''p'' の値による $100p$ パーセンタイルを求める関数である。
ここでは ''p'' は ''0.3'' としているので、30 パーセンタイルをしきい値としている。

----
==== k-means法 ====
データ分析の方法である非階層的クラスター分析の __**k-means法**__ で画素値を二つのクラスターに自動的に分類して、しきい値を自動で決める方法である。

この方法では、ヒストグラムを二つのクラスターに適当に分け、各画素値をクラスターの平均に距離的に近い方のクラスターに入れるということを繰り返してクラスターを確定させる。 
そして、そのクラスターの境界をしきい値とする。

{{ detection_k-means.png?nolink |k-means法}}

k-means 法は、以下で述べる ''threshold(im,thr)'' で ''thr'' を ''%%"auto"%%'' とすればできる。

===== 2値化の関数 =====

imager で2値化を行う関数として ''threshold()'' がある。
<code rsplus>
> threshold(cimg変数, thr)
</code>
ただし、cimg 変数はグレースケール画像でなくてはならない。

''thr'' には以下のしきい値を指定する。

^  しきい値        ^  機能                                                                   ^
| ''%%"auto"%%''   | k-means 法によって自動的にしきい値が決められる（デフォルト）。          |
| 実数             | 実数の場合は、モード法におけるしきい値になる。                          |
| ''%%"100p%"%%''  | p-タイル法における $100p$ % になる。''100p'' のところには数字を入れる。  |

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

次のプログラム k-means.R は画像を k-means法で2値化する。

<file rsplus k-means.R[enable_line_numbers="true",highlight_lines_extra="9"]>
# k-means法

library(imager)

# グレースケール画像
im <- grayscale(load.image("https://www.kitp.org/numericalipb/nameplate.jpg"))

# k-means法による2値化画像
im.bin <- threshold(im, "auto")

layout(t(1 : 2))
plot(im, main = "元の画像")
plot(im.bin, main = "2値化")
layout(1)
</file>

9 行目の ''%%"auto"%%'' を ''0.5'' のように実数にすればモード法、''%%"50%"%%'' のようにパーセントを文字列で入れれば p-タイル法による2値化になる。

===== 2値画像の平滑化 =====

threshold.R でしきい値 ''thr'' を ''0.3'' にして実行してみると、白い部分に小さな黒い点もしくは黒い文字の上に小さな白い点が見えるだろう。
これは 2値画像でよく現れるノイズで、その見た目から__**ごま塩ノイズ**__という。

{{ nameplate_binary.png?nolink&600 |2値画像}}

以前に行った平滑化でこのノイズを除去できるが、2値画像では 2値であることを利用して__**膨張・収縮**__ または __**収縮・膨張**__ によってノイズを除去する方法がある。

imager で膨張と収縮はそれぞれ以下の意味になる。
  * 膨張：白の画素を膨張させること。
  * 収縮：黒の画素を膨張させること（結果として白の画素が収縮する）。

imager 以外では膨張と収縮が逆（黒と白とが逆）の意味の場合があるので注意してほしい。

imager で膨張と収縮を行う関数は以下の通り。

^  関数                              ^  機能                                                           ^
| ''dilate_square(cimg変数, size)''  | 画素を正方形に膨張する。 ''size'' は縦横に膨張するピクセル数である。  |
| ''erode_square(cimg変数, size)''   | 画素を正方形に収縮する。''size'' は縦横に収縮するピクセル数である。   |


----
==== 膨張・収縮 ====
膨張したあとに収縮することで<fc #ff0000>**黒のノイズを除去**</fc>する。

=== ①膨張 ===
まず白の画素を膨張する。
以下の図では白の画素一つ一つを 3x3 の正方形に膨張している。
これによって白い部分が膨張して黒い孤立したノイズは除去される。

=== ②収縮 ===
その後、黒の画素を膨張する。
以下の図では黒の画素一つ一つを 3x3 の正方形に膨張している。
結果として白い部分が収縮して残った部分が元のサイズに戻る。

{{ dilate_erode.png?nolink |膨張・収縮の例}}

次のプログラム binary_de.R は膨張・収縮を行う。
膨張と収縮のサイズは 3 にしている。

<file rsplus binary_de.R[enable_line_numbers="true",highlight_lines_extra="18,21"]>
# 膨張・収縮

library(imager)

# グレースケール画像
im <- grayscale(load.image("https://www.kitp.org/numericalipb/nameplate.jpg"))

# しきい値
thr <- 0.3

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

# imの画素値がしきい値以上なら画素値1とする
im.bin[im >= thr] <- 1

# 膨張
im.bin1 <- dilate_square(im.bin, 3)

# 収縮
im.bin2 <- erode_square(im.bin1, 3) 

layout(t(matrix(1 : 4, 2, 2)))
plot(im, interpolate = FALSE, main = "元の画像")
plot(im.bin, interpolate = FALSE, main = "2値化")
plot(im.bin1, interpolate = FALSE, main = "膨張")
plot(im.bin2, interpolate = FALSE, main = "収縮")
layout(1)
</file>

----
==== 収縮・膨張 ====
収縮したあとに膨張することで<fc #ff0000>**白のノイズを除去**</fc>する。

=== ①収縮 ===
まず黒の画素を膨張する。
以下の図では黒の画素一つ一つを 3x3 の正方形に膨張している。
結果として白い部分が収縮して白い孤立したノイズは除去される。

=== ②膨張 ===
その後、白の画素を膨張する。
以下の図では白の画素一つ一つを 3x3 の正方形に膨張している。
白い部分が膨張して残った部分が元のサイズに戻る。


{{ erode_dilate.png?nolink |収縮・膨張}}

binary_de.R で、''delate_square()'' と ''erode_square()'' を入れ替えれば、収縮・膨張になる。

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

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


===== 物体検出とラベリング =====

画像から必要な物体だけを取り出したい場合は、まずその物体がどこからどこまでかを検出しなければならない。
画像から物体を検出するためには、濃淡や色の違い、形の違いなどで検出することが考えられる。
最近はAIの機械学習による物体検出が流行であるが、ここでは取り扱わない。

ここでは2値画像を使って、画素値が 1 （白）のひとまとまりのつながった領域（__**連結領域**__）を認識する方法で物体検出を行う。
連結領域に同じ番号でラベル付けをし、異なる図形（物体）を分類する処理を__**ラベリング**__という。

ラベリングをした後に、ラベルを指定すればそのラベルの付いた図形のみを抜き出すことができる。

ラベリングは関数 ''label()'' で行う。
<code rsplus>
> label(cimg変数)
</code>
それぞれの図形に番号のラベル付けが行われ、その図形に属する画素の画素値をその番号にした cimg 変数が返ってくる。
例えば、番号 5 が付けられた図形に含まれる画素の画素値はすべて 5 になる。

次のプログラム label.R はサンプル画像の ''coins'' を k-means 法で2値化し、ラベリングを行って図形を抜き出すプログラムである。

<file rsplus label.R[enable_line_numbers="true"]>
# ラベリング

library(imager)

# グレースケール画像
im <- grayscale(load.example("coins"))

# 2値化
im.bin <- threshold(im, "auto")

# ラベリング
x <- label(im.bin)

# ラベル付けされた図形の面積（画素数）
areas <- table(x)

# 面積が100以上の図形のラベルと面積を出力
print(areas[areas >= 100])

# ラベルを指定して図形を抜き出す
im.ex <- im
im.ex[x != 113 & x != 238] <- 1

layout(t(1 : 2))
plot(im, main = "元画像")
plot(im.ex, main = "抜き出した図形")
layout(1)
</file>

12行目で 2 値画像がラベリングされ、返ってきた cimg 変数を ''x'' に代入している。

15行目で各ラベルの画素数（面積）をクロス集計している。

18行目で面積が 100 以上の図形のラベルと面積を出力している。
出力は以下の通りである。
黄色の行がラベルで、その下の行が面積である。
<code rsplus [highlight_lines_extra="2,4"]>
x
    1     4    66    92    96    98   113   118   231   233   235   236   238 
 8755 69597  2459  1684   123  1631  1193  1133  1834  1325  1203  1133  1129 
  244   292   300   302   305   311   313   425   426   434   439   461   475 
 1104  3054  1633  1352  1461  1095  1148  2099  1954  1918  1728  1312  1462 
</code>

このプログラムでは以下のように各コインにラベル付けされる（ラベル 1 のコインは背景との分離が不十分）。
{{ label_coin.png?nolink |ラベルの番号}}

20行目以降でラベルを指定して画像を抜き出している。
22行目ではラベル「113」「238」の図形以外の画素値を 1（白）にしている。


===== 課題 =====

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

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