【画像処理】NumpyでLoGフィルターとゼロ交差によるエッジ検出

NumpyでLoG(Laplacian of Gaussian)フィルターを実装して、ゼロ交差によるエッジ検出を試してみます。

LoGフィルターとは名前の通りガウス関数のラプラシアン(二次微分)で、以下の式で表されます。

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

まず、使用する画像を読み込んで、カラー画像をグレースケール画像に変換します。

```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
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/82be8cb0-2e8f-4e0d-b7ab-c540d0f43788.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)
```

LoGフィルターのカーネルを作成して、先ほどのグレースケール画像に畳み込みます。

```python
def create_log_kernel(size=(13, 13), sigma=2):
  center = ((size[0] - 1) / 2, (size[1] - 1) / 2)
  sigma2 = sigma * sigma
  sigma6 = sigma2 * sigma2 * sigma2
  sigma2_2 = sigma2 * 2.0
  sigma2_2_inv = 1.0 / sigma2_2
  sigma6_2pi_inv = 1.0 / (sigma6 * 2.0 * np.pi)
  def calc_weight(y, x):
    sqDist = (x - center[1]) ** 2 + (y - center[0]) ** 2
    return (sqDist - sigma2_2) * sigma6_2pi_inv * np.exp(-sqDist * sigma2_2_inv)
  return np.fromfunction(calc_weight, size)

log_kernel = create_log_kernel()
log_image = convolve2d(gray_image, log_kernel)
value_range = max(abs(log_image.min()), abs(log_image.max()))
plt.imshow(log_image, cmap='bwr', vmin=-value_range, vmax=value_range)
plt.colorbar()
```

![LogImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/572fb4fc-9b06-789b-6be0-0c6188266f21.png)

LoGフィルターを畳み込んだ画像からゼロ交差を求めて、エッジを検出します。ここではある画素について左右または上下の画素の正負が異なるときに、その画素がゼロ交差とみなしています。

```python
def detect_zero_crossing(image):
  pad_image = np.pad(image, 1, 'edge')
  shape = image.shape + (3, 3)
  strides = pad_image.strides * 2
  strided_image = np.lib.stride_tricks.as_strided(pad_image, shape, strides).reshape(shape[0], shape[1], 9)
  return np.apply_along_axis(
      lambda array: 1 if array[3] * array[5] < 0 or array[2] * array[6] < 0 else 0,
      2, strided_image)

zero_crossing_image = detect_zero_crossing(log_image)
plt.imshow(zero_crossing_image, cmap='gray')
```

![ZeroCrossing.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/8f049a37-784a-7549-25e2-0cfd59731bcc.png)

ゼロ交差までできたので、いくつかパラメータを変更して試してみます。基本的には、`sigma`が大きくなるほど、`size`も大きくする必要があります。

```python
log_kernel1 = create_log_kernel(size=(7, 7), sigma=1)
log_image1 = convolve2d(gray_image, log_kernel1)
value_range = max(abs(log_image1.min()), abs(log_image1.max()))
plt.imshow(log_image1, cmap='bwr', vmin=-value_range, vmax=value_range)
plt.colorbar()
```

![LogImage1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/33173d00-ae02-0c1a-deb2-9bf6d8efa4a2.png)

```python
zero_crossing_image1 = detect_zero_crossing(log_image1)
plt.imshow(zero_crossing_image1, cmap='gray')
```

![ZeroCrossing1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/f80ae56d-3e10-cf50-b2c5-d102e12ac328.png)

```python
log_kernel2 = create_log_kernel(size=(31, 31), sigma=4)
log_image2 = convolve2d(gray_image, log_kernel2)
value_range = max(abs(log_image2.min()), abs(log_image2.max()))
plt.imshow(log_image2, cmap='bwr', vmin=-value_range, vmax=value_range)
plt.colorbar()
```

![LogImage2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/49b1a83a-15d3-2d7a-2198-d4fc8fbeba06.png)

```python
zero_crossing_image2 = detect_zero_crossing(log_image2)
plt.imshow(zero_crossing_image2, cmap='gray')
```

![ZeroCrossing2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/5da2a9e8-fb7d-1c23-58e6-b5285477b060.png)


```python
log_kernel3 = create_log_kernel(size=(41, 41), sigma=7)
log_image3 = convolve2d(gray_image, log_kernel3)
value_range = max(abs(log_image3.min()), abs(log_image3.max()))
plt.imshow(log_image3, cmap='bwr', vmin=-value_range, vmax=value_range)
plt.colorbar()
```

![LogImage3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/6b4100c6-e689-1398-c31d-6640d3a0727c.png)

```python
zero_crossing_image3 = detect_zero_crossing(log_image3)
plt.imshow(zero_crossing_image3, cmap='gray')
```

![ZeroCrossing3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/1be105c6-bc73-f6fe-0199-c52e9bfebddc.png)

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

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