# 画像処理 エッジ検出器(1)

前回の空間フィルタリングに続いて、エッジ検出についての理論とPythonによる実装方法を掲載しています。
コードに間違いなどがあればご指摘いただけると幸いです。

#### エッジ検出器とは
エッジ検出器とは、画像に写っている物体の輪郭を検出する処理のことです。

身近な例で言うと、コロナ禍になってZoomなどのオンラインでの会話が増えましたが、背景以下の写真のように変化させる手法も、エッジ検出の応用例です。

![zoom.jpeg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1726209/1dba20a9-884c-55ad-ded0-f6ed01882d0b.jpeg)
写真引用元 : [デザイン情報サイト「JDN」](https://www.japandesign.ne.jp/news/2020/04/53234/)

本記事では、このようにモノの輪郭を検出する基礎的なエッジ検出器について記していきます。

#### エッジ検出の適用例

ここでは、白黒画像の明るさが急激な変化をしているところを物体の縁として検出することを目的とします。

例として、一次元で明るさが以下のグラフの$(a)$のように変化している場合を考えます。
グラフの$(a)$のグラフについて、一次微分したグラフを$(b)$、二次微分したグラフを$(c)$とします。
![graph.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1726209/1f6639ee-0644-fc87-63fe-ca39ca54d97e.png)
すると、$(b)$のグラフでは、色が急激に変化している箇所のみが高い値となること、そして$(c)$のグラフでは、同じところにおいて符号が変化するゼロ交差点があることがわかります。

このように、一次微分の計算または二次微分の計算をすることによって、エッジ検出が可能であることが考えられます。

本記事では主に一次微分のエッジ検出について記しております。二次微分を利用したエッジ検出では次回以降の記事にで記していこうと思います。


#### 一次微分を利用したエッジ検出の理論

位置$x$における$f(x)$の微分の式は以下のように表されます。

```math
f'(x) = \lim_{\Delta x \to 0} \frac{f(x + \Delta x) - f(x)}{\Delta x}
```

連続値の場合の微分の式はこのようになりますが、実際に画素について扱う場合は離散的なグラフとなりますので、$\Delta x = ±1$として、$f(x ± 1)$を位置$x$に隣り合う画素値とすると、以下のように単純な差分についての式が得られます。

```math
\begin{align}
f'(x) &= f(x + 1) - f(x) \\
f'(x) &= f(x) - f(x - 1)
\end{align}
```
この二式から$f(x)$を消去すると、以下の式が得られます。

```math
f'(x) = \frac{f(x + 1) - f(x - 1)}{2}
```

上記のことを$I(x, y)$について適用すると、$x$方向と$y$方向の微分の式$I_x$と$I_y$はそれぞれ以下のように表されます。
```math
\begin{align}
I_x &= \frac{I(x + 1, y) - I(x - 1, y)}{2} \\
I_y &= \frac{I(x, y + 1) - I(x, y - 1)}{2}
\end{align}
```
この式から、注目画素$I(x, y)$に対して、隣り合った画素の差分を$2$で割ったものを計算し、それを可視化することによってエッジ検出が可能であることがわかります。
![calculate.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1726209/379ad4da-c427-1fb6-868e-4b976b2ecd76.png)
上記の式は以下のカーネルを適用することと同値です。
```math
F_x = \frac{1}{2}
\begin{pmatrix}
0 & 0 & 0 \\
-1 & 0 & 1 \\
0 & 0 & 0 
\end{pmatrix}
,
\hspace{10px}
F_y = \frac{1}{2}
\begin{pmatrix}
0 & -1 & 0 \\
0 & 0 & 0 \\
0 & 1 & 0 
\end{pmatrix}
```
$x$方向と$y$方向の検出結果をまとめるには、各画素の微分値の二乗平均平方根を計算します。
```math
∇I = \sqrt{I_x^2 + I_y^2}
```
ただ、このままでは薄い縁まで検出されてしまうため、この$∇I$に閾値を設けて信頼性の高い縁のみを残す処理(Hysteresis Threshold処理)必要になります。また、上記の画像のようにエッジが太くなってしまうためそれを細くする処理(Non maximum Suppression処理)が必要になります。後述するPythonではこれらの処理は施しておりませんが、後ほど記事にしていく予定です。

また、上記のカーネルではノイズに弱いため、注目画素の周辺の画素値の微分も計算し、それらの値の平均を取る方法が考えられています。代表的なものがプリューウィットフィルタ(Prewwit Filter)、ソーベルフィルタ(Sobel Filter)が挙げられます。

##### プリューウィットフィルタ(Prewwit Filter)
```math
F_x =
\begin{pmatrix}
-1 & 0 & 1 \\
-1 & 0 & 1 \\
-1 & 0 & 1 
\end{pmatrix}
,
\hspace{10px}
F_y =
\begin{pmatrix}
-1 & -1 & -1 \\
0 & 0 & 0 \\
1 & 1 & 1 
\end{pmatrix}
```

##### ソーベルフィルタ(Sobel Filter)
```math
F_x =
\begin{pmatrix}
-1 & 0 & 1 \\
-2 & 0 & 2 \\
-1 & 0 & 1 
\end{pmatrix}
,
\hspace{10px}
F_y =
\begin{pmatrix}
-1 & -2 & -1 \\
0 & 0 & 0 \\
1 & 2 & 1 
\end{pmatrix}
```

#### キャニーエッジ検出器(Canny Edge)とは
キャニーエッジ検出器とは、前回の記事で説明したガウシアンフィルタと上記のソーベルフィルタを組み合わせたエッジ検出器です。
上記のプリューウィットフィルタとソーベルフィルタはノイズ対策としてカーネルを再検討された結果ですが、それでもノイズの影響を避けるには限界があります。
それを防ぐために、ノイズ対策として、前回の記事で説明したガウシアンフィルタを施した後にソーベルフィルタでエッジ検出する方法が考えられます。
これによってノイズに強い適切なエッジ検出が期待できます。

#### エッジ検出のPythonでの実装

##### 1.プリューウィットフィルタ(Prewwit Filter) & ソーベルフィルタ(Sobel Filter)

この二つはカーネルを変更するだけでできますので、同じソースコードで実装します。

##### Pythonでの実装

```python
import cv2
import numpy as np

# フィルタ計算する関数
def filter2D(src, kernel):
    # カーネルサイズ
    m, n = kernel.shape
    # 畳み込み演算をしない領域の幅
    ignore = int((m - 1) / 2)
    height, width = src.shape
    
    img_edge = np.zeros((height, width))
    
    for y in range(ignore, height - ignore):
        for x in range(ignore, width - ignore):
            # 畳み込み演算
            img_edge[y][x] = np.sum(src[y - ignore: y + ignore + 1, x - ignore: x + ignore + 1] * kernel)
    
    return img_edge

# 画像を白黒画像としてインポート
img_gray = cv2.imread('rankuru_noise.jpg', cv2.IMREAD_GRAYSCALE)

# プリューウィットフィルタのカーネルの設定
prewwit_kernel_x = np.array([
    [-1,  0,  1],
    [-1,  0,  1],
    [-1,  0,  1],
])
prewwit_kernel_y = np.array([
    [-1, -1, -1],
    [ 0,  0,  0],
    [ 1,  1,  1],
])

# ソーベルフィルタのカーネルの設定
sobel_kernel_x = np.array([
    [-1,  0,  1],
    [-1,  0,  1],
    [-1,  0,  1],
])
sobel_kernel_y = np.array([
    [-1, -1, -1],
    [ 0,  0,  0],
    [ 1,  1,  1],
])


# プリューウィットフィルタを適用
prewwit_img_gray_x = filter2D(img_gray, prewwit_kernel_x)
prewwit_img_gray_y = filter2D(img_gray, prewwit_kernel_y)
# ソーベルフィルタを適用
sobel_img_gray_x = filter2D(img_gray, sobel_kernel_x)
sobel_img_gray_y = filter2D(img_gray, sobel_kernel_y)

# x,yの検出結果をまとめる
prewwit_img_edge = np.sqrt(prewwit_img_gray_x ** 2, prewwit_img_gray_x ** 2)
sobel_img_edge = np.sqrt(sobel_img_gray_x ** 2, sobel_img_gray_x ** 2)
```

```filter_2D```関数を使用することで、カーネルに基づいた画素の計算ができるようになっています。
また、各画素の微分値の二乗平均平方根を計算してまとめます。
これによってソーベルフィルタとプリューウィットフィルタを適用させることができます。

計算結果を表示してみます。```matplotlib.pyplot```モジュールの```imshow```関数を使用します。

```python
import matplotlib.pyplot as plt

# 元の白黒画像を表示
plt.figure(figsize=(7, 6))
plt.imshow(img_gray, cmap='gray')
plt.xticks([]);plt.yticks([])
plt.show()

# フィルタ適用後の画像の表示
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(14, 6))
ax = ax.flatten()
ax[0].imshow(prewwit_img_edge, cmap='gray_r')
ax[1].imshow(sobel_img_edge, cmap='gray_r')
ax[0].set_xticks([]);ax[0].set_yticks([]);ax[0].set_title('Prewwit Filter')
ax[1].set_xticks([]);ax[1].set_yticks([]);ax[1].set_title('Sobel Filter')
plt.show()
```
得られた画像は以下のようになりました。
* 元の画像
![gray.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1726209/7a719bf0-2898-da52-a9ed-0c11a252db42.png)
* フィルタ適用後の画像
![ps.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1726209/bdde3cc8-6ea3-070c-33a2-043a4422e610.png)
プリューウィットフィルタとソーベルフィルタの双方に大きな差はない結果となりました。
また、双方ともランドクルーザの大まかな輪郭がかたどれていることがわかります。

##### 2.キャニーエッジ検出器(Canny Edge)

キャニーエッジ検出器も、上記のソーベルフィルタと前回の記事のガウシアンフィルタを組み合わせるだけなので、特に難しくはないかと思います。

##### Pythonでの実装

```python
import cv2
import numpy as np

# キャニーエッジ検出器

# 1. ガウシアンフィルタ

# カーネル関数
def kernel_function(x, y, sigma):
    z = (1 / (2 * np.pi * sigma ** 2)) * np.exp(-(x ** 2 + y ** 2) / (2 * sigma ** 2))
    return z

# ガウシアンフィルタ用のカーネルを作成する関数
def create_kernel(size=(3, 3), sigma=1):
    kernel = np.zeros(size)
    kerner_x, kernel_y = int((size[1] - 1) / 2), int((size[0] - 1) / 2)
    for x in range(size[1]):
        for y in range(size[0]):
            kernel[y][x] = kernel_function(x - kerner_x, y - kernel_y, sigma)
    # カーネルの規格化
    kernel_norm = kernel / np.sum(kernel)
    return kernel_norm

# ガウシアンフィルタの計算
def gaussian_filter(img, kernel):
    # カーネルサイズを取得
    line, column = kernel.shape
    # フィルタをかける画像の高さと幅を取得
    height, width = img.shape

    # 畳み込み演算をしない領域の幅を指定
    height_ignore = int((line - 1) / 2)
    width_ignore = int((column - 1) / 2)

    # 出力画像用の配列
    img_filter = img.copy()

    # フィルタリングの計算
    for y in range(height_ignore, height - height_ignore):
        for x in range(width_ignore, width - width_ignore):
            img_filter[y][x] = np.sum(img[y - height_ignore: y + height_ignore + 1, x - width_ignore: x + width_ignore + 1] * kernel)
    return img_filter

# 2. ソーベルフィルタ

# ソーベルフィルタの計算
def sobel_filter(src, kernel):
    # カーネルサイズ
    line, column = kernel.shape
    # 畳み込み演算をしない領域の幅
    ignore = int((line - 1) / 2)
    height, width = src.shape
    
    img_edge = np.zeros((height, width))
    
    for y in range(ignore, height - ignore):
        for x in range(ignore, width - ignore):
            # 畳み込み演算
            img_edge[y][x] = np.sum(src[y - ignore: y + ignore + 1, x - ignore: x + ignore + 1] * kernel)
    
    return img_edge

# 入力画像を読み込み
img_gray = cv2.imread("rankuru_noise.jpg", cv2.IMREAD_GRAYSCALE)

# カーネルを定義(行数と列数は奇数である必要あり)
gaussian_kernel = create_kernel(size=(7, 7), sigma=5)

# ガウシアンフィルタを適用
img_filter = gaussian_filter(img_gray, gaussian_kernel)

# ソーベルフィルタのカーネルを定義
sobel_kernel_x = np.array([
    [-1,  0,  1],
    [-2,  0,  2],
    [-1,  0,  1],
])
sobel_kernel_y = np.array([
    [-1, -2, -1],
    [ 0,  0,  0],
    [ 1,  2,  1],
])

# ソーベルフィルタを適用
img_edge_x = sobel_filter(img_filter, sobel_kernel_x)
img_edge_y = sobel_filter(img_filter, sobel_kernel_y)

# x,y方向に分けて計算した結果を一つに集約
img_edge = np.sqrt(img_edge_x ** 2, img_edge_y ** 2)
```
上記のコードでキャニーエッジ検出器が適用できます。
計算結果を表示してみます。
```python
import matplotlib.pyplot as plt
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(14, 6))
ax = ax.flatten()
ax[0].imshow(img_gray, cmap='gray')
ax[1].imshow(img_edge, cmap='gray_r')
ax[0].set_xticks([]);ax[0].set_yticks([])
ax[1].set_xticks([]);ax[1].set_yticks([])
plt.show()
```
元画像と得られた画像の比較は以下のようになります。
* キャニーエッジ検出器
![canny.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1726209/e1fb7d1f-597b-e9be-2c60-12fa6d1b5730.png)
フィルタリング処理を施しているため、他のものよりもノイズが少ない結果となりました。

#### まとめ

エッジ検出についてと、その代表例である「プリューウィットフィルタ」と「ソーベルフィルタ」、さらにフィルタリングと組み合わせた「キャニーエッジ検出器」のPythonによる自力実装を行いました。
次回は二次微分を用いたエッジ検出について記していこうと思います。

以上になります。