【画像処理】Numpyで平滑化フィルター

Numpyで画像の平滑化フィルターを実装してみます。

まず、画像の読み込みを行います。

```python
original_image = plt.imread(image_name)
if np.issubdtype(original_image.dtype, np.integer):
  original_image = original_image / np.iinfo(original_image.dtype).max
plt.imshow(original_image)
```

![OriginalImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/7ec8e4d5-a9ea-8ee7-d1ac-c298c1d8cad4.png)

画像の畳み込みを行う関数を定義しておきます。詳細は以下の記事を参照してください。
[【画像処理】Numpyで空間フィルタリング\(畳み込み演算\) \- Qiita](https://qiita.com/aa_debdeb/items/e74eceb13ad8427b16c6)

```python
def _convolve2d(image, kernel):
  shape = (image.shape[0] - kernel.shape[0] + 1, image.shape[1] - kernel.shape[1] + 1) + kernel.shape
  strides = image.strides * 2
  strided_image = np.lib.stride_tricks.as_strided(image, shape, strides)
  return np.einsum('kl,ijkl->ij', kernel, strided_image)

def _convolve2d_multichannel(image, kernel):
  convolved_image = np.empty((image.shape[0] - kernel.shape[0] + 1, image.shape[1] - kernel.shape[1] + 1, image.shape[2]))
  for i in range(image.shape[2]):
    convolved_image[:,:,i] = _convolve2d(image[:,:,i], kernel)
  return convolved_image

def _pad_singlechannel_image(image, kernel_shape, boundary):
  return np.pad(image, ((int(kernel_shape[0] / 2),), (int(kernel_shape[1] / 2),)), boundary)

def _pad_multichannel_image(image, kernel_shape, boundary):
  return  np.pad(image, ((int(kernel_shape[0] / 2),), (int(kernel_shape[1] / 2),), (0,)), boundary)

def convolve2d(image, kernel, boundary='edge'):
  if image.ndim == 2:
    pad_image = _pad_singlechannel_image(image, kernel.shape, boundary) if boundary is not None else image
    return _convolve2d(pad_image, kernel)
  elif image.ndim == 3:
    pad_image = _pad_multichannel_image(image, kernel.shape, boundary) if boundary is not None else image
    return _convolve2d_multichannel(pad_image, kernel)
```

## 平均化フィルター

近傍画素での平均をとるようにカーネル内の重みが同じで合計が1になるフィルターを平均化フィルターと呼びます。
例えば、3x3の場合および5x5の場合はそれぞれ次のようになります。

```math
\begin{bmatrix}
\frac{1}{9} & \frac{1}{9} & \frac{1}{9} \\
\frac{1}{9} & \frac{1}{9} & \frac{1}{9} \\
\frac{1}{9} & \frac{1}{9} & \frac{1}{9} \\
\end{bmatrix}
```

```math
\begin{bmatrix}
\frac{1}{25} & \frac{1}{25} & \frac{1}{25} & \frac{1}{25} & \frac{1}{25} \\
\frac{1}{25} & \frac{1}{25} & \frac{1}{25} & \frac{1}{25} & \frac{1}{25} \\
\frac{1}{25} & \frac{1}{25} & \frac{1}{25} & \frac{1}{25} & \frac{1}{25} \\
\frac{1}{25} & \frac{1}{25} & \frac{1}{25} & \frac{1}{25} & \frac{1}{25} \\
\frac{1}{25} & \frac{1}{25} & \frac{1}{25} & \frac{1}{25} & \frac{1}{25} \\
\end{bmatrix}
```

平均化フィルターのカーネルを作成する関数は次のようになります。

```python
def create_averaging_kernel(size = (5, 5)):
  return np.full(size, 1 / (size[0] * size[1]))
```

パラメータを変えて平均化フィルター適用してみます。

```python
averaging_kernel_3x3 = create_averaging_kernel(size=(3, 3))
averaging_image_3x3 = convolve2d(original_image, averaging_kernel_3x3)
plt.imshow(averaging_image_3x3)
```

![AveragingFilter3x3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/cf2be2ff-a3ac-69cd-beb7-fe3bce96004e.png)

```python
averaging_kernel_5x5 = create_averaging_kernel()
averaging_image_5x5 = convolve2d(original_image, averaging_kernel_5x5)
plt.imshow(averaging_image_5x5)
```

![AveragingFilter5x5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/1e2335fa-cebc-baa1-3536-6831d6b5f350.png)

```python
averaging_kernel_11x11 = create_averaging_kernel(size=(11, 11))
averaging_image_11x11 = convolve2d(original_image, averaging_kernel_11x11)
plt.imshow(averaging_image_11x11)
```

![AveragingFilter11x11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/60a16291-25c5-d191-ca2c-6b6987a3c9ae.png)

縦方向だけで平均化してみます。

```python
averaging_kernel_17x1 = create_averaging_kernel(size=(17, 1))
averaging_image_17x1 = convolve2d(original_image, averaging_kernel_17x1)
plt.imshow(averaging_image_17x1)
```

![AveragingFilter17x1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/89557e71-13fa-5adc-a6aa-0c4e6f7aa670.png)


## ガウシアンフィルター

中心画素に近いほど影響が大きくなるように、以下のガウス関数に基づき重みづけをしたフィルターをガウシアンフィルターと呼びます。

```math
f(x, y) = \frac{1}{2\pi\sigma^2}exp(-\frac{x^2+y^2}{2\sigma^2})
```


ガウシアンフィルターのカーネルを作成する関数は次のようになります。重みの合計が1になるように最後に正規化しています。

```python
def create_gaussian_kernel(size=(5, 5), sigma=1):
  center = ((size[0] - 1) / 2, (size[1] - 1) / 2)
  sigma2 = 2 * sigma * sigma
  kernel = np.fromfunction(lambda y, x: np.exp(-((x - center[1]) ** 2 + (y - center[0]) ** 2) / sigma2), size)
  kernel = kernel / np.sum(kernel)
  return kernel
```

パラメータを変えてガウシアンフィルターを適用してみます。

```python
gaussian_kernel1 = create_gaussian_kernel()
gaussian_image1 = convolve2d(original_image, gaussian_kernel1)
plt.imshow(gaussian_image1)
```

![GaussianFilter1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/3891e5c5-76c3-e260-2eb2-d4f7b147d8be.png)

```python
gaussian_kernel2 = create_gaussian_kernel(size=(9, 9), sigma=3)
gaussian_image2 = convolve2d(original_image, gaussian_kernel2)
plt.imshow(gaussian_image2)
```

![GaussianFilter2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/41b6b26f-7f8a-071a-1788-d3c5c1e1fe6a.png)

```python
gaussian_kernel3 = create_gaussian_kernel(size=(15, 15), sigma=5)
gaussian_image3 = convolve2d(original_image, gaussian_kernel3)
plt.imshow(gaussian_image3)
```

![GaussianFilter3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/2f72aab2-679a-77bc-4378-74f13fa2ba92.png)

---

実装したコードはGoogle Colaboratoryに置いておきました。

https://colab.research.google.com/drive/1WiUTYobuJrSWopnhwC7iMhChz9OLLUDY?usp=sharing