【画像処理】NumpyでCannyエッジ検出

NumpyでCannyエッジ検出を実装してみます。

## 画像読み込み

画像を読み込み、カラー画像をグレースケール画像に変換します。

```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/9220c8e3-4253-4c2a-ec60-644d5997cbff.png)

## 畳み込み関数

畳み込みを行う関数`convolve2d`を定義しておきます。詳細は以下の記事を参照してください。
[【画像処理】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)
```

## 画像の平滑化(ノイズ除去)

ガウシアンフィルターで画像を平滑化して、ノイズを除去します。

```python
def create_gaussian_kernel(size=(3, 3), 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 

gaussian_kernel = create_gaussian_kernel()
gaussian_image = convolve2d(gray_image, gaussian_kernel)
plt.imshow(gaussian_image, cmap='gray')
```

![GaussianImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/5ddce9e0-d7c1-7535-0695-ac8a8c8b8a82.png)

## 画像の微分

Sobelフィルターで画像の微分を行い、勾配を求めます。

まず横方向の微分です。

```python
sobel_x_kernel = np.array([[-1, 0, 1],
                           [-2, 0, 2],
                           [-1, 0, 1]])
sobel_x_image = convolve2d(gaussian_image, sobel_x_kernel) / 8
value_range = max(abs(sobel_x_image.min()), abs(sobel_x_image.max()))
plt.imshow(sobel_x_image, cmap='bwr', vmin=-value_range, vmax=value_range)
plt.colorbar()
```

![SobelXImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/35aa3923-741e-d2d0-9b24-db1fc65ac16e.png)

次に縦方向の微分です。

```python
sobel_y_kernel = np.array(([[1, 2, 1],
                            [0, 0, 0],
                            [-1, -2, -1]]))
sobel_y_image = convolve2d(gaussian_image, sobel_y_kernel) / 8
value_range = max(abs(sobel_y_image.min()), abs(sobel_y_image.max()))
plt.imshow(sobel_y_image, cmap='bwr', vmin=-value_range, vmax=value_range)
plt.colorbar()
```

![SobelYImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/70172d5b-312b-7c76-4d0c-6d0e46248342.png)

各方向の微分から勾配の強度と方向を求めます。

```python
sobel_intensity_image = np.sqrt(sobel_x_image ** 2 + sobel_y_image ** 2)
plt.imshow(sobel_intensity_image, cmap='gray')
plt.colorbar()
```

![SobelIntensityImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/a580a47a-813c-db3a-9498-5d7fa95ce690.png)

```python
sobel_direction_image = np.arctan(np.divide(sobel_y_image, sobel_x_image, out=np.zeros_like(sobel_y_image), where=sobel_x_image!=0))
plt.imshow(sobel_direction_image, cmap='hsv')
plt.colorbar()
```

![SobelDirectionImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/e22d4ee2-da19-772b-1ceb-90ab64503dc7.png)

## 勾配方向の非極値抑制


細線化のために、先ほど求めた勾配強度のうち勾配方向で最大値となるものだけを残します。

```python
def _suppress_nonmaximum(array):
  direction = array[9]
  if direction < 0:
    direction += np.pi
  if direction < 0.125 * np.pi:
    v0, v1, v2 = (array[3], array[4], array[5])
  elif direction < 0.375 * np.pi:
    v0, v1, v2 = (array[2], array[4], array[6])
  elif direction < 0.625 * np.pi:
    v0, v1, v2 = (array[1], array[4], array[7])
  elif direction < 0.875 * np.pi:
    v0, v1, v2 = (array[0], array[4], array[8])
  else:
    v0, v1, v2 = (array[3], array[4], array[5])
  return v1 if v1 >= v0 and v1 >= v2 else 0

def suppress_nonmaximum(intensity_image, direction_image):
  pad_image = np.pad(intensity_image, 1, 'edge')
  shape = intensity_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(_suppress_nonmaximum, 2, np.dstack((strided_image, direction_image)))

nonmaximum_suppression_image = suppress_nonmaximum(sobel_intensity_image, sobel_direction_image)
plt.imshow(nonmaximum_suppression_image, cmap='gray')
plt.colorbar()
```

![SuppressNonmaximumImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/3e8314be-c5d9-8dba-62d5-c463e3a7a644.png)

## ヒステリシス閾値処理

非極値抑制した画像からエッジ画像を作成します。まず、値が閾値$T_{high}$よりも大きい画素をエッジとみなします。さらに値が閾値$T_{low}$よりも大きく、かつエッジが隣接している画素もエッジとみなすことで連結したエッジを検出します。

```python
def histeresis_thresholding(image, threshold):
  edge_image = np.where(image >= threshold[1], 1.0, 0.0)

  while True:
    pad_edge_image = np.pad(edge_image, 1, 'edge')
    shape = image.shape + (3, 3)
    strides = pad_edge_image.strides * 2
    neighbour_edge_image = np.lib.stride_tricks.as_strided(pad_edge_image, shape, strides)
    edge_candidate_index = np.where((image >= threshold[0]) & (edge_image != 1.0))

    changed = False
    for index in list(zip(*edge_candidate_index)):
      if np.any(neighbour_edge_image[index] == 1.0):
        edge_image[index] = 1.0
        changed = True
    if not changed:
      break;

  return edge_image

canny_image = histeresis_thresholding(nonmaximum_suppression_image, (0.03, 0.1))
plt.imshow(canny_image, cmap='gray')
```

![HisteresisThresholdingImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/e35b8122-49c9-7a34-7d76-ba499ce80e59.png)

## まとめ

以上の処理をまとめて、Cannyエッジ検出を行う関数を作成すると次のようになります。

```python
def canny(image, gaussian_size=(3, 3), gaussian_sigma=1, threshold=(0.02, 0.05)):
  # Noise reduction by gaussian filter
  gaussian_kernel = create_gaussian_kernel(size=gaussian_size, sigma=gaussian_sigma)
  gaussian_image = convolve2d(gray_image, gaussian_kernel)

  # Gradient by sobel filter
  sobel_x_kernel = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
  sobel_x_image = convolve2d(gaussian_image, sobel_x_kernel) / 8
  sobel_y_kernel = np.array(([[1, 2, 1], [0, 0, 0], [-1, -2, -1]]))
  sobel_y_image = convolve2d(gaussian_image, sobel_y_kernel) / 8
  sobel_intensity_image = np.sqrt(sobel_x_image ** 2 + sobel_y_image ** 2)
  sobel_direction_image = np.arctan(np.divide(sobel_y_image, sobel_x_image, out=np.zeros_like(sobel_y_image), where=sobel_x_image!=0))

  # Non-maximum suppression
  nonmaximum_suppression_image = suppress_nonmaximum(sobel_intensity_image, sobel_direction_image)

  # Histresis thresholding
  return histeresis_thresholding(nonmaximum_suppression_image, threshold)
```

パラメータを変えてCannyエッジ検出を試してみます。

```python
canny_image1 = canny(gray_image)
plt.imshow(canny_image1, cmap='gray')
```

![CannyImage1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/45d5ea2d-9ca1-4320-af1f-dfadb8436471.png)

```python
canny_image2 = canny(gray_image, gaussian_size=(5, 5), gaussian_sigma=0.2, threshold=(0.01, 0.04))
plt.imshow(canny_image2, cmap='gray')
```

![CannyImage2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/600eb486-86f3-38fe-7316-86ce9f76790c.png)

```python
canny_image3 = canny(gray_image, gaussian_size=(3, 3), gaussian_sigma=1.2, threshold=(0.06, 0.08))
plt.imshow(canny_image3, cmap='gray')
```

![CannyImage3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/113b69b7-c195-d8b3-208f-89bf94bc0231.png)

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

https://colab.research.google.com/drive/1zg9dcdljN5VT3N6gAVwa4XiHMYJ-_LX3?usp=sharing