【画像処理】Numpyでテンプレートマッチング

Numpyでテンプレートマッチングを実装してみます。

まず、マッチングを行う画像を読み込みます。

```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/27bc81e7-869e-af9f-8b31-1b98db273230.png)

テンプレート画像を用意します。

```python
template_image = gray_image[10:55,60:120]
plt.figure(figsize=(1, 1))
plt.imshow(template_image, cmap='gray')
```

![TemplateImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/f0cb83fe-0099-1ebf-9afd-3ad131bad561.png)

## SSD

SSD(Sum of Squared Difference)で相違度を計算します。

```math
R_{ssd}(x, y) = \sum_{i=0}^{N-1}\sum_{j=0}^{M-1}(I(x+i,y+j) - T(i,j))^2
```

$R_{ssd}(x,y)$は位置$(x,y)$におけるSSDの値、$I(i,j)$と$T(i,j)$はそれぞれ位置$(i,j)$におけるマッチング対象の画像の画素値とテンプレートの画素値、$(N,M)$はテンプレートの大きさとなります。

```python
def match_template_ssd(image, template):
  shape = (image.shape[0] - template.shape[0] + 1, image.shape[1] - template.shape[1] + 1) + template.shape
  strided_image = np.lib.stride_tricks.as_strided(image, shape, image.strides * 2)
  return np.sum((strided_image - template) ** 2.0, axis=(2, 3))

ssd_image = match_template_ssd(gray_image, template_image)
plt.imshow(ssd_image, cmap='gray')
plt.colorbar()
```

![SsdImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/228c6af4-679a-0bb2-0abd-f5ae0477ed72.png)

SSDが最も小さい箇所がテンプレートと最もマッチする箇所になります。そこに四角形を描画しています。四角形の描画に関しては以下の記事を参照してください。
[【画像処理】Numpyで図形描画 \- Qiita](https://qiita.com/aa_debdeb/items/8178fef706e8dc7e30d2)

```python
def stroke_rectangle(image, color, weight, center, size):
  coord = np.fromfunction(lambda y, x: np.dstack((y + 0.5, x + 0.5)), image.shape[:2])
  dist = np.dstack((abs(coord[:,:,1] - center[0]) - size[0] / 2, abs(coord[:,:,0] - center[1]) - size[1] / 2)).max(axis=2)
  condition = abs(dist) - weight * 0.5 <= 0
  if image.ndim == 3:
    condition = np.tile(condition.reshape(condition.shape + (1,)), (1, 1, image.shape[2]))
  return np.where(condition, color, image)

index = np.unravel_index(np.argmin(ssd_image), ssd_image.shape)
ssd_rect_image = stroke_rectangle(original_image, np.array([1.0, 0.0, 1.0]), 2.0, (index[1] + template_image.shape[1] / 2, index[0] + template_image.shape[0] / 2), (template_image.shape[1], template_image.shape[0]))
plt.imshow(ssd_rect_image)
```

![SsdMatchImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/bc09edc2-3d3b-f6fd-b594-ec9d3e02283a.png)

## SAD

SAD(Sum of Absolute Difference)で相違度を計算します。

```math
R_{sad}(x, y) = \sum_{i=0}^{N-1}\sum_{j=0}^{M-1}|I(x+i,y+j) - T(i,j)|
```

```python
def match_template_sad(image, template):
  shape = (image.shape[0] - template.shape[0] + 1, image.shape[1] - template.shape[1] + 1) + template.shape
  strided_image = np.lib.stride_tricks.as_strided(image, shape, image.strides * 2)
  return np.sum(np.abs(strided_image - template), axis=(2, 3))

sad_image = match_template_sad(gray_image, template_image)
plt.imshow(sad_image, cmap='gray')
plt.colorbar()
```

![SadImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/949e4576-f89d-2578-deaf-fc48cc323b32.png)

SSDと同様にSADの値が最も小さい箇所が最もマッチする箇所なので、そこに四角形を描画します。

```python
index = np.unravel_index(np.argmin(sad_image), sad_image.shape)
sad_rect_image = stroke_rectangle(original_image, np.array([1.0, 0.0, 1.0]), 2.0, (index[1] + template_image.shape[1] / 2, index[0] + template_image.shape[0] / 2), (template_image.shape[1], template_image.shape[0]))
plt.imshow(sad_rect_image)
```

![SadMatchImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/c0ae30b0-66af-5071-41cb-d93ccdd78133.png)

## NCC

NCC(Normalized Cross-Correlation)で類似度を計算します。

```math
R_{ncc}(x, y) = \frac{\sum_{i=0}^{N-1}\sum_{j=0}^{M-1}I(x+i,y+j)T(i,j)}{\sqrt{\sum_{i=0}^{N-1}\sum_{j=0}^{M-1}I(x+i,y+j)^2\times\sum_{i=0}^{N-1}\sum_{j=0}^{M-1}T(i,j)^2}}
```

```python
def match_template_ncc(image, template):
  shape = (image.shape[0] - template.shape[0] + 1, image.shape[1] - template.shape[1] + 1) + template.shape
  strided_image = np.lib.stride_tricks.as_strided(image, shape, image.strides * 2)
  return np.sum(strided_image * template, axis=(2, 3)) \
    / (np.sqrt(np.sum(strided_image * strided_image, axis=(2, 3)) * np.sum(template * template)))

ncc_image = match_template_ncc(gray_image, template_image)
plt.imshow(ncc_image, cmap='gray')
plt.colorbar()
```

![NccImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/f8e9bd2f-f774-ca6b-0009-fdad9aa5d068.png)

NCCは最も値が大きい箇所が最もマッチする箇所になるので、そこに四角形を描画します。

```python
index = np.unravel_index(np.argmax(ncc_image), ncc_image.shape)
ncc_rect_image = stroke_rectangle(original_image, np.array([1.0, 0.0, 1.0]), 2.0, (index[1] + template_image.shape[1] / 2, index[0] + template_image.shape[0] / 2), (template_image.shape[1], template_image.shape[0]))
plt.imshow(ncc_rect_image)
```

![NccMatchImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/31491101-a096-9c7a-2d00-10496036f4ec.png)

---

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

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