【画像処理】Numpyでバイラテラルフィルター

Numpyでバイラテラル(Bilateral)フィルターを実装してみます。

ガウシアンフィルターでは空間方向の距離に応じてガウス分布で加重平均を求めることで画像を平滑化しました。バイラテラルフィルターでは空間方向の距離に加えて、画素値の差に応じてもガウス分布で重み付けします。これにより、エッジを保存したまま平滑化によりノイズを除去することができます。

入力画像を$f(i, j)$、出力画像を$g(i, j)$とすると、次のように表現できます。

```math
g(i,j)=\frac{\sum_{n}\sum_{m}w(i,j,m,n)f(i+m,j+n)}{\sum_{n}\sum_{m}w(i,j,m,n)} \\
w(i,j,m,n)=exp(-\frac{m^2+n^2}{2\sigma^2_{space}})exp(-\frac{(f(i,j)-f(i+m,j+n))^2}{2\sigma^2_{pixel}})
```

---

まず、検証用のテスト画像を作成します。

```python
import numpy as np
import matplotlib.pyplot as plt

test_image = np.fromfunction(lambda y, _: np.where(y < 128, 1.0, 0.0), (256, 256))
test_image = np.clip(test_image + np.random.normal(0, 0.2, test_image.shape), 0.0, 1.0)
plt.imshow(test_image, cmap='gray')
```

![TestImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/5523fd04-633c-11a5-5609-0fffa00713d3.png)

Bilateralフィルターを適用する関数を定義します。

```python
def bilateral_filter(image, size=(5, 5), space_sigma=1, pixel_sigma=0.2, boundary='edge'):
  pad_image = np.pad(image, ((int(size[0] / 2),), (int(size[1] / 2),)), boundary)
  areas = np.lib.stride_tricks.as_strided(pad_image, image.shape + size, pad_image.strides * 2)
  centers = np.tile(image.reshape(image.shape + (1, 1)), (1, 1) + size)
  dists = np.fromfunction(lambda y, x: (x - int(size[1] / 2)) ** 2 + (y - int(size[0] / 2)) ** 2, size)
  weights = np.exp(-dists / (2.0 * space_sigma * space_sigma)) * np.exp(-((centers - areas) ** 2) / (2.0 * pixel_sigma * pixel_sigma * 2.0))
  weight_sum = np.sum(weights, axis=(2, 3))
  return np.einsum('ijkl,ijkl->ij', weights, areas) / weight_sum
```

バイラテラルフィルターをテスト画像に繰り返し適用します。

1回目

```python
bilateral_image1 = bilateral_filter(test_image)
plt.imshow(bilateral_image1, cmap='gray')
```

![BilateralImage1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/4d9d3e70-7890-2065-1997-e8245ef2de91.png)

2回目

```python
bilateral_image2 = bilateral_filter(bilateral_image1)
plt.imshow(bilateral_image2, cmap='gray')
```

![BilateralImage2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/8b47f02f-3628-0956-7356-b65430156756.png)

3回目

```python
bilateral_image3 = bilateral_filter(bilateral_image2)
plt.imshow(bilateral_image3, cmap='gray')
```

![BilateralImage3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/cebfc305-9233-fa44-b685-ba328f79f901.png)

エッジを維持したまま、ノイズが低減することを確認できました。

---

比較のためにガウシアンフィルターをテスト画像に適用してみます。ガウシアンフィルターについては以下の記事を参照してください。
[【画像処理】Numpyで平滑化フィルター \- Qiita](https://qiita.com/aa_debdeb/items/9c5d5f93890265abfdec)

```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)

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[0]) ** 2 + (y - center[1]) ** 2) / sigma2), size)
  kernel = kernel / np.sum(kernel)
  return kernel

gaussian_kernel = create_gaussian_kernel()
```

1回目

```python
gaussian_image1 = convolve2d(test_image, gaussian_kernel)
plt.imshow(gaussian_image1, cmap='gray')
```

![GaussianImage1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/7358bf05-c2c3-3fc6-5c8d-338dc7615b12.png)

2回目

```python
gaussian_image2 = convolve2d(gaussian_image1, gaussian_kernel)
plt.imshow(gaussian_image2, cmap='gray')
```

![GaussianImage2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/801aa3ef-d19b-189f-ed0d-30cab941a3a5.png)

3回目

```python
gaussian_image3 = convolve2d(gaussian_image2, gaussian_kernel)
plt.imshow(gaussian_image3, cmap='gray')
```

![GaussianImage3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/73e02661-9045-86c7-a25d-9d741f4a9236.png)

バイラテラルフィルターの結果と比較して、エッジが鈍っていることが確認できます。

---

バイラテラルフィルターをカラー画像に適用してみます。
バイラテラルフィルターをカラー画像に繰り返し適用すると、イラスト調の画像になります。

使用する画像を読み込みます。

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

![ColorImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/a9f24d6e-0a44-4b13-b512-4a100dfb1790.png)

マルチチャネル画像にバイラテラルフィルターを適用する関数を定義します。

```python
def bilateral_filter_multichannel(image, size=(5, 5), space_sigma=1, pixel_sigma=0.2, boundary='edge'):
  bilateral_image = np.empty(image.shape)
  for i in range(image.shape[2]):
    bilateral_image[:,:,i] = bilateral_filter(image[:,:,i], size, space_sigma, pixel_sigma, boundary)
  return bilateral_image 
```

バイラテラルフィルターを繰り返し適用します。

1回目

```python
bilateral_color_image1 = bilateral_filter_multichannel(color_image, size=(11, 11), space_sigma=3, pixel_sigma=0.05)
plt.imshow(bilateral_color_image1)
```

![ColorBilateralImage1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/4c2c2073-91c3-331e-0724-476a5b35ca00.png)

2回目

```python
bilateral_color_image2 = bilateral_filter_multichannel(bilateral_color_image1, size=(11, 11), space_sigma=3, pixel_sigma=0.05)
plt.imshow(bilateral_color_image2)
```

![ColorBilateralImage2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/2af20ad0-f5eb-f79f-f995-2b3bbd711488.png)

3回目

```python
bilateral_color_image3 = bilateral_filter_multichannel(bilateral_color_image2, size=(11, 11), space_sigma=3, pixel_sigma=0.05)
plt.imshow(bilateral_color_image3)
```

![ColorBilateralImage3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/e89d45ce-654b-99cb-35dc-281c2da8ad8b.png)


---

実装したコードはGoogle Colaboratoryに置いてあります。

https://colab.research.google.com/drive/1gFHYDRvuvy_hQGcSHKhyYq6IZa4Kf-jK?usp=sharing