【画像処理】NumpyでHarisのコーナー検出

NumpyでHarisのコーナー検出を実装してみます。

コーナーを検出する画像を読み込みます。

```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/c644edc1-81df-76f3-738b-e3efe842a9a7.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)
```

Sobelフィルターで横方向、縦方向の勾配画像$I_x$、$I_y$をそれぞれ作成します。

```python
sobel_x_kernel = np.array([[-1, 0, 1],
                           [-2, 0, 2],
                           [-1, 0, 1]])
sobel_x_image = convolve2d(gray_image, sobel_x_kernel)
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/433bd0d0-9a5e-b90a-9f6c-e773dac5166e.png)

```python
sobel_y_kernel = np.array(([[1, 2, 1],
                            [0, 0, 0],
                            [-1, -2, -1]]))
sobel_y_image = convolve2d(gray_image, sobel_y_kernel)
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/e09c8686-c6b7-78aa-5411-aa8a8ea66752.png)


勾配画像から$I_{xx}=I_x\cdot I_x$、$I_{yy}=I_y\cdot I_y$、$I_{xy}=I_x\cdot I_y$を求めます。

```python
sobel_x2_image = sobel_x_image ** 2
plt.imshow(sobel_x2_image, cmap='gray')
plt.colorbar()
```

![SobelX2Image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/ab2edb49-e927-9f13-0d09-d161274a11f0.png)

```python
sobel_y2_image = sobel_y_image ** 2
plt.imshow(sobel_y2_image, cmap='gray')
plt.colorbar()
```

![SobelY2Image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/cdf528f6-acf5-f09d-8247-7a73d7bf110d.png)

```python
sobel_xy_image = sobel_x_image * sobel_y_image
value_range = max(abs(sobel_xy_image.min()), abs(sobel_xy_image.max()))
plt.imshow(sobel_xy_image, cmap='bwr', vmin=-value_range, vmax=value_range)
plt.colorbar()
```

![SobelXyImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/b2fdd7b0-7cf0-cc4d-428b-8193033eb085.png)

$I_{xx}$、$I_{yy}$、$I_{xy}$から局所領域で平均化した画像$S_{xx}$、$S_{yy}$、$S_{xy}$をそれぞれ作成します。ここではガウシアンフィルターを使用しています。ガウシアンフィルターについては以下の記事を参照してください。
[【画像処理】Numpyで平滑化フィルター \- Qiita](https://qiita.com/aa_debdeb/items/9c5d5f93890265abfdec)

```python
gaussian_sobel_x2_image = convolve2d(sobel_x2_image, gaussian_kernel)
plt.imshow(gaussian_sobel_x2_image, cmap='gray')
plt.colorbar()
```

![GaussianSobelX2Image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/984c6f86-c5e6-e259-6b0d-80da283d9142.png)

```python
gaussian_sobel_y2_image = convolve2d(sobel_y2_image, gaussian_kernel)
plt.imshow(gaussian_sobel_y2_image, cmap='gray')
plt.colorbar()
```

![GaussianSobelY2Image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/0df0e94b-a80c-77c7-7708-95339d875d18.png)

```python
gaussian_sobel_xy_image = convolve2d(sobel_xy_image, gaussian_kernel)
value_range = max(abs(gaussian_sobel_xy_image.min()), abs(gaussian_sobel_xy_image.max()))
plt.imshow(gaussian_sobel_xy_image, cmap='bwr', vmin=-value_range, vmax=value_range)
plt.colorbar()
```

![GaussianSobelXYImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/f0dcc5dd-d4ae-6072-a686-8cd3c07d8932.png)

画素$(x,y)$における行列$\boldsymbol{M}$を次のように定義して$R$を求めます。この$R$が一定値以上であれば、その画素をコーナーとみなすことができます。

```math
\boldsymbol{M}(x,y)=
\begin{bmatrix}
S_{xx}(x,y) & S_{xy}(x,y) \\
S_{xy}(x,y) & S_{yy}(x,y) \\
\end{bmatrix} \\
R=det\boldsymbol{M} - k(tr\boldsymbol{M})^2
```

```python
m = np.stack([np.dstack((gaussian_sobel_x2_image, gaussian_sobel_xy_image)),
              np.dstack((gaussian_sobel_xy_image, gaussian_sobel_y2_image))], axis=-1)
det_m = np.linalg.det(m)
tr_m = np.trace(m, axis1=2, axis2=3)
k = 0.06
r = det_m - k * (tr_m ** 2)
plt.imshow(r)
plt.colorbar()
```

![RImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/8f7c5370-6924-1312-921f-624294abcf62.png)

$R$が一定以上の値の箇所を抽出して、そこに円を描画します。円の描画については以下の記事を参照してください。
[【画像処理】Numpyで図形描画 \- Qiita](https://qiita.com/aa_debdeb/items/8178fef706e8dc7e30d2)

```python
def fill_circle(image, color, center, radius):
  coord = np.fromfunction(lambda y, x: np.dstack((y + 0.5, x + 0.5)), image.shape[:2])
  dist2 = (coord[:,:,1] - center[0]) ** 2 + (coord[:,:,0] - center[1]) ** 2
  condition = dist2 <= radius ** 2
  if image.ndim == 3:
    condition = np.tile(condition.reshape(condition.shape + (1,)), (1, 1, image.shape[2]))
  return np.where(condition, color, image)

corners = np.where(r > 1.0)
corner_image = original_image
for coord in list(zip(corners[1], corners[0])):
  corner_image = fill_circle(corner_image, np.array([1.0, 0.0, 1.0]), coord, 3.0)
plt.imshow(corner_image)
```

![CornerImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/25ed4ca0-2316-906d-bff3-7901bfde045b.png)

1つのコーナー付近で複数の検出があるので、$R$が局所最大となる画素のみを残します。

```python
def is_local_maximum(image, size=(5, 5), 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, r.shape + size, pad_image.strides * 2)
  local_maximum = np.max(areas, axis=(2, 3))
  return image == local_maximum

corners = np.where(is_local_maximum(r) & (r > 1.0))
corner_image = original_image
for coord in list(zip(corners[1], corners[0])):
  corner_image = fill_circle(corner_image, np.array([1.0, 0.0, 1.0]), coord, 3.0)
plt.imshow(corner_image)
```

![CornerLocalMaxImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/d0b27357-1dab-a608-78e8-072b9c4908bf.png)

以上、Harisのコーナー検出によりコーナーを検出することができました。

---

これまでの操作をまとめて一つの関数にします。

```python
def detect_corners_by_haris(image, threshold, k=0.06, gaussian_size=(5, 5), gaussian_sigma=1, local_maximum_size=(5, 5)):
  dx = convolve2d(image, np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]]))
  dy = convolve2d(image, np.array(([[1, 2, 1], [0, 0, 0], [-1, -2, -1]])))
  dx2 = dx ** 2
  dy2 = dy ** 2
  dxy = dx * dy
  gauss_kernel = create_gaussian_kernel(size=gaussian_size, sigma=gaussian_sigma)
  gdx2 = convolve2d(dx2, gauss_kernel)
  gdy2 = convolve2d(dy2, gauss_kernel)
  gdxy = convolve2d(dxy, gauss_kernel)
  m = np.stack([np.dstack((gdx2, gdxy)), np.dstack((gdxy, gdy2))], axis=-1)
  det_m = np.linalg.det(m)
  tr_m = np.trace(m, axis1=2, axis2=3)
  r = det_m - k * (tr_m ** 2)
  return np.where(is_local_maximum(r, size=local_maximum_size) & (r > threshold))

def draw_circles(image, coordinates, color, radius):
  for coordinate in list(zip(coordinates[1], coordinates[0])):
    image = fill_circle(image, color, coordinate, radius)
  return image

corners = detect_corners_by_haris(gray_image, 1.0)
corner_image = draw_circles(original_image, corners, np.array([1.0, 0.0, 1.0]), 3.0)
plt.imshow(corner_image)
```

![CornerLocalMaxImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/d0b27357-1dab-a608-78e8-072b9c4908bf.png)

---

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

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