【画像処理】NumpyでDoGフィルター

NumpyでDoG(Diffrence of Gaussian)フィルターを実装してみます。
DoGフィルターでは、$\sigma$の値が異なる2つのガウシアンフィルターをそれぞれ適用し、その差分を求めます。

まず、DoGフィルターを適用する画像を読み込んでグレースケール画像にします。

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

original_image = plt.imread(image_name)
if np.issubdtype(original_image.dtype, np.integer):
  original_image = original_image / np.iinfo(original_image.dtype).max
gray_image = 0.2116 * original_image[:,:,0] + 0.7152 * original_image[:,:,1] + 0.0722 * original_image[:,:,2]
plt.imshow(gray_image, cmap='gray')
```

![GrayImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/cec845c4-39f5-3f10-c85c-f65b720ae2c0.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[1]) ** 2 + (y - center[0]) ** 2) / sigma2), size)
  kernel = kernel / np.sum(kernel)
  return kernel
```

DoGフィルターを適用します。

```python
def dog_filter(image, sigma1=2.0, sigma2=1.0, size=(7, 7), boundary='edge'):
  kernel1 = create_gaussian_kernel(size=size, sigma=sigma1)
  kernel2 = create_gaussian_kernel(size=size, sigma=sigma2)
  gauss_image1 = convolve2d(image, kernel1, boundary=boundary)
  gauss_image2 = convolve2d(image, kernel2, boundary=boundary)
  return gauss_image1 - gauss_image2

dog_image = dog_filter(gray_image)
value_range = max(abs(dog_image.min()), abs(dog_image.max()))
plt.imshow(dog_image, cmap='bwr', vmin=-value_range, vmax=value_range)
plt.colorbar()
```

![DogImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/88b840e9-c4e8-f022-b18f-10979d982c4c.png)

---

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

https://colab.research.google.com/drive/1h4ep-9-SBiJN-C8tyaCirUeEFSm6-fDn?usp=sharing