python+opencvで画像処理の勉強5 幾何学的変換

pythonとopencvを使って画像処理を勉強していきます。

前回
python+opencvで画像処理の勉強4 周波数領域におけるフィルタリング
周波数フィルタリングについて学びました。

https://qiita.com/tanaka_benkyo/items/bfa35e7f08faa7b7a985

今回は幾何学的変換について学び、最後は応用として画像のつなぎ合わせを行います。

使用画像は旭川動物園のペンギンの散歩の写真です。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/622cdc12-faba-c532-6c83-6818cde69498.png)

## 線形変換
### 線形変換の一般系
座標$(x,y)$の位置の点が、変換により座標$(x',y')$に移動するとする。
そのとき、以下の式で表される変換を**線形変換**と呼ぶ。

```math
    \left[\begin{array}{c}
    x'\\
    y'
    \end{array}\right] =
    \left[\begin{array}{cc}
    a&b\\
    c&d
    \end{array}\right]
    \left[\begin{array}{c}
    x\\
    y
    \end{array}\right]
```

### 拡大・縮小
**拡大・縮小**は、以下の式で表される。

```math
    \left[\begin{array}{c}
    x'\\
    y'
    \end{array}\right] =
    \left[\begin{array}{cc}
    s_{x}&0\\
    0&s_{y}
    \end{array}\right]
    \left[\begin{array}{c}
    x\\
    y
    \end{array}\right]
```

ここで$s_{x},s_{y}$はx方向、y方向の拡大（縮小）率である。

```python
width = img_rgb.shape[1]
height = img_rgb.shape[0]

x_scale = [2,1,3]
y_scale = [2,2,1]

fig, ax = plt.subplots(1, 3, figsize=(16, 9))
for i in range(3):
    dst = cv2.resize(img_rgb, (width*x_scale[i], height*y_scale[i]))
    ax[i].imshow(dst)
    ax[i].set_title('x_scale={} y_scale={}'.format(str(x_scale[i]), str(y_scale[i])))
    ax[i].set_xticks([])
    ax[i].set_yticks([])
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/424f377f-8245-7471-c5ec-e374fe68222a.png)


### 回転
原点を中心に、反時計回りに角度$\theta$だけ**回転**する変換は、以下の式で表される。

```math
    \left[\begin{array}{c}
    x'\\
    y'
    \end{array}\right] =
    \left[\begin{array}{cc}
    \cos{\theta}&-\sin{\theta}\\
    \sin{\theta}&\cos{\theta}
    \end{array}\right]
    \left[\begin{array}{c}
    x\\
    y
    \end{array}\right]
```

ここでは、画像の左下を中心に回転した例を示す。

```python
fig, ax = plt.subplots(1, 2, figsize=(10, 5))

center = (0, height)

affine_trans = cv2.getRotationMatrix2D(center, 20.0, 1.0)
dst = cv2.warpAffine(img_rgb, affine_trans, (width, height))
ax[0].imshow(dst);
ax[0].set_xticks([]);
ax[0].set_yticks([]);

affine_trans = cv2.getRotationMatrix2D(center, -30.0, 1.0)
dst = cv2.warpAffine(img_rgb, affine_trans, (width, height), flags=cv2.INTER_CUBIC)
ax[1].imshow(dst);
ax[1].set_xticks([]);
ax[1].set_yticks([]);
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/f294e66a-94c9-ee3a-493d-f2f890cef5d1.png)


### 鏡映
ある直線に対して、対象な位置に反転する変換を**鏡映**と呼ぶ。
x軸に関する変換

```math
    \left[\begin{array}{c}
    x'\\
    y'
    \end{array}\right] =
    \left[\begin{array}{cc}
    1&0\\
    0&-1
    \end{array}\right]
    \left[\begin{array}{c}
    x\\
    y
    \end{array}\right]
```

y軸に関する変換

```math
    \left[\begin{array}{c}
    x'\\
    y'
    \end{array}\right] =
    \left[\begin{array}{cc}
    -1&0\\
    0&1
    \end{array}\right]
    \left[\begin{array}{c}
    x\\
    y
    \end{array}\right]
```

直線y=xに関する変換

```math
    \left[\begin{array}{c}
    x'\\
    y'
    \end{array}\right] =
    \left[\begin{array}{cc}
    0&1\\
    1&0
    \end{array}\right]
    \left[\begin{array}{c}
    x\\
    y
    \end{array}\right]
```

```python
fig, ax = plt.subplots(1, 4, figsize=(20, 5))

ax[0].imshow(img_rgb)
ax[0].set_xticks([])
ax[0].set_yticks([])

flip_n = [0, 1, -1]

for i in range(1,4):
    dst = cv2.flip(img_rgb, flip_n[i-1])
    ax[i].imshow(dst)
    ax[i].set_title('flip={}'.format(str(flip_n[i-1])))
    ax[i].set_xticks([])
    ax[i].set_yticks([])
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/c7660a6c-b6d8-adbc-626b-09ce271b1fbb.png)


### スキュー
長方形を傾けて平行四辺形にするような変換を**スキュー**または**せん断**と呼ぶ。
x軸方向へのスキューは以下の式で表される。

```math
    \left[\begin{array}{c}
    x'\\
    y'
    \end{array}\right] =
    \left[\begin{array}{cc}
    1&b\\
    0&1
    \end{array}\right]
    \left[\begin{array}{c}
    x\\
    y
    \end{array}\right]
```

傾ける角度を$\theta$とすると、$b=\tan{\theta}$の関係がある。
また、y軸方向へのスキューは以下の式で表される。

```math
    \left[\begin{array}{c}
    x'\\
    y'
    \end{array}\right] =
    \left[\begin{array}{cc}
    1&0\\
    c&1
    \end{array}\right]
    \left[\begin{array}{c}
    x\\
    y
    \end{array}\right]
```

傾ける角度を$\theta$とすると、$c=\tan{\theta}$の関係がある。

```python
fig, ax = plt.subplots(1, 2, figsize=(10, 5))

mat = np.array([[1, 0, 0], [0.3, 1, 0]], dtype=np.float32)
dst = cv2.warpAffine(img_rgb, mat, (width, int(height + width * 0.3)))
ax[0].imshow(dst);
ax[0].set_xticks([]);
ax[0].set_yticks([]);

mat = np.array([[1, 0.6, 0], [0, 1, 0]], dtype=np.float32)
dst = cv2.warpAffine(img_rgb, mat, (int(width + height * 0.6), height))
ax[1].imshow(dst);
ax[1].set_xticks([]);
ax[1].set_yticks([]);
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/722feeb0-5da3-0d8b-1163-e38f14a1e0f8.png)


### 合成変換
ここまでに示した各種変換は、ベクトルと行列を記号で置き換えることで、以下のように表される。

```math
    \boldsymbol{x'}=\boldsymbol{A}\boldsymbol{x}
```

ここで、

```math
    \boldsymbol{x}=
    \left[\begin{array}{c}
    x\\
    y
    \end{array}\right],
    \boldsymbol{x'}=
    \left[\begin{array}{c}
    x'\\
    y'
    \end{array}\right],
    \boldsymbol{A}=
    \left[\begin{array}{cc}
    a&b\\
    c&d
    \end{array}\right]
```

線形変換は2×2の行列$\boldsymbol{A}$で表現される。
座標$x$が変換$\boldsymbol{A}$により、座標$\boldsymbol{x'}$に移ったとすると以下のように表される。

```math
    \boldsymbol{x'}=\boldsymbol{A}\boldsymbol{x}
```

さらに、座標$\boldsymbol{x'}$が、変換$\boldsymbol{B}$により座標$\boldsymbol{x''}$に移ったとすると、以下のように表される。

```math   \boldsymbol{x''}=\boldsymbol{B}\boldsymbol{x'}
```

以下の式が導かれる。

```math
    \boldsymbol{x''}=(\boldsymbol{BA})\boldsymbol{x}=\boldsymbol{Cx}
```

ここで、$\boldsymbol{C}$は$\boldsymbol{C}=\boldsymbol{BA}$なる2×2行列であり、これも1つの線形変換を表すことになる。
任意の線形変換の組み合わせは、やはり線形変換となる。一般に、個々の順番を入れ替えると、結果が異なることに注意する。

## 同次座標とアフィン変換・射影変換
### 平行移動
x軸方向の移動量$t_{x}$、y軸方向の移動量を$t_{y}$としたとき、この**平行移動**は以下の式で表される。

```math
    x'=x+t_{x}\\
    y'=y+t_{y}
```

単純なものであるが、線形変換の一般系で表現することはできない。

```python
fig, ax = plt.subplots(1, 2, figsize=(10, 5))

mat = np.array([[1, 0, 200], [0, 1, 50]], dtype=np.float32)
dst = cv2.warpAffine(img_rgb, mat, (width, height))
ax[0].imshow(dst);
ax[0].set_xticks([]);
ax[0].set_yticks([]);

mat = np.array([[1, 0, -100], [0, 1, -50]], dtype=np.float32)
dst = cv2.warpAffine(img_rgb, mat, (width, height))
ax[1].imshow(dst);
ax[1].set_xticks([]);
ax[1].set_yticks([]);
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/9b199fc5-604f-9b85-3df3-c732e98cf481.png)


### 同次座標
座標$(x,y)$に対し、要素を1つ増やした座標$(\xi_{1},\xi_{2},\xi_{3})$を、以下の関係式を満たすように定義する。

```math
    x=\frac{\xi_{1}}{\xi_{3}}\\
    y=\frac{\xi_{2}}{\xi_{3}}
```

ただし、$\xi_{1},\xi_{2},\xi_{3}$の少なくとも1つは0でないとする。
このように定義される座標を**同次座標**と呼ぶ。
同次座標において、$\lambda\neq0$となる任意の$\lambda$に対し、$(\lambda\xi_{1},\xi_{2},\lambda\xi_{3})$と$(\xi_{1},\lambda\xi_{2},\xi_{3})$は,
通常の座標に直したとき、ともに$(\frac{\xi_{1}}{\xi_{3}},\frac{\xi_{2}}{\xi_{3}})$となる。
同次座標系においては、定数倍をしても変わらないとみなせる。
このような関係を**同値**であると呼び、以下のように表す。

```math
    \left[\begin{array}{c}
    \xi_{1}\\
    \xi_{2}\\
    \xi_{3}
    \end{array}\right] \sim
    \left[\begin{array}{c}
    \lambda\xi_{1}\\
    \lambda\xi_{2}\\
    \lambda\xi_{3}
    \end{array}\right]
```

ベクトル表記では、頭に~を付けて表す。
$$\boldsymbol{\tilde{\xi}}\sim\lambda\boldsymbol{\tilde{\xi}}$$

### アフィン変換
同次座標を利用すると、平行移動は以下の式のように表すことができる。

```math
    \left[\begin{array}{c}
    x'\\
    y'\\
    1
    \end{array}\right] \sim
    \left[\begin{array}{ccc}
    1&0&t_{x}\\
    0&1&t_{y}\\
    0&0&1
    \end{array}\right]
    \left[\begin{array}{c}
    x\\
    y\\
    1
    \end{array}\right]
```

また、線形変換は、同次座標を用いると以下のように表現される。

```math
    \left[\begin{array}{c}
    x'\\
    y'\\
    1
    \end{array}\right] \sim
    \left[\begin{array}{ccc}
    a&b&0\\
    c&d&0\\
    0&0&1
    \end{array}\right]
    \left[\begin{array}{c}
    x\\
    y\\
    1
    \end{array}\right]
```

さらに、任意の線形変換と平行移動を組み合わせた変換を**アフィン変換**と呼ぶ。
アフィン変換の一般系は以下のように求めることができる。

```math
    \left[\begin{array}{ccc}
    1&0&t_{x}\\
    0&1&t_{y}\\
    0&0&1
    \end{array}\right]
    \left[\begin{array}{ccc}
    a&b&0\\
    c&d&0\\
    0&0&1
    \end{array}\right] \sim
    \left[\begin{array}{ccc}
    a&b&t_{x}\\
    c&d&t_{y}\\
    0&0&1
    \end{array}\right]
```

```python
a = [1, 2, 3]
b = [4, -2, -5]
c = [4, 0.5, 2]
d = [-3, 2, 0]

tx = [10, 1400, 1000]
ty = [30, 0, 200]

fig, ax = plt.subplots(1, 3, figsize=(16, 5))

for i in range(3):
    mat = np.array([[a[i], b[i], tx[i]], [c[i], d[i], ty[i]]], dtype=np.float32)
    dst = cv2.warpAffine(img_rgb, mat, (7*width, 7*height))
    ax[i].imshow(dst)
    ax[i].set_xticks([])
    ax[i].set_yticks([])
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/ebafea38-36ad-3aa0-93f4-45da0bf9132f.png)

アフィン変換の特殊な例で、任意の回転と平行移動を組みわせたものを**ユークリッド変換**と呼ぶ。
対象の大きさを保ったまま、任意に移動するような変換である

```math
    \left[\begin{array}{ccc}
    \cos{\theta}&-\sin{\theta}&t_{x}\\
    \sin{\theta}&\cos{\theta}&t_{y}\\
    0&0&1
    \end{array}\right]
```

```python
thetas = [np.pi/6, -np.pi/3, np.pi/4]
c = np.cos(thetas)
s = np.sin(thetas)

tx = [20, 70, 50]
ty = [-100, 150, -150]

fig, ax = plt.subplots(1, 3, figsize=(16, 5))

for i in range(3):
    mat = np.array([[c[i], -s[i], tx[i]], [s[i], c[i], ty[i]]], dtype=np.float32)
    dst = cv2.warpAffine(img_rgb, mat, (width, height))
    ax[i].imshow(dst)
    ax[i].set_xticks([])
    ax[i].set_yticks([])
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/b0bf7004-d7ef-a510-023f-99e24b0bac3b.png)

ユークリッド変換に対し、縦横の倍率が等しい拡大・縮小を加えたものを**相似変換**と呼ぶ。

```math
    \left[\begin{array}{ccc}
    s\cos{\theta}&-s\sin{\theta}&t_{x}\\
    s\sin{\theta}&s\cos{\theta}&t_{y}\\
    0&0&1
    \end{array}\right]
```

```python
thetas = [np.pi/6, -np.pi/3, np.pi/4]
c = np.cos(thetas)
s = np.sin(thetas)
p = [4, 2, 0.6]

tx = [20, 150, 50]
ty = [-150, 100, -150]

fig, ax = plt.subplots(1, 3, figsize=(16, 5))

for i in range(3):
    mat = np.array([[p[i]*c[i], -p[i]*s[i], tx[i]], [p[i]*s[i], p[i]*c[i], ty[i]]], dtype=np.float32)
    dst = cv2.warpAffine(img_rgb, mat, (width, height))
    ax[i].imshow(dst)
    ax[i].set_xticks([])
    ax[i].set_yticks([])
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/e4353e04-f480-9307-e78b-f7cbf53e68f6.png)


### 射影変換
同次座標を用いて、さらに一般的な変換を以下の式で表現することもできる。
この変換を**射影変換（ホモグラフィ）**と呼ぶ。

```math
    \left[\begin{array}{c}
    x'\\
    y'\\
    1
    \end{array}\right] \sim
    \left[\begin{array}{ccc}
    h_{11}&h_{12}&h_{13}\\
    h_{21}&h_{22}&h_{23}\\
    h_{31}&h_{32}&h_{33}
    \end{array}\right]
    \left[\begin{array}{c}
    x\\
    y\\
    1
    \end{array}\right]
```

これをベクトルと行列の記号を用いて表すと、
$$\boldsymbol{\tilde{x}'}\sim\boldsymbol{H}\boldsymbol{\tilde{x}}$$
ただし、

```math
    \boldsymbol{\tilde{x}}=
    \left[\begin{array}{c}
    x\\
    y\\
    1
    \end{array}\right],
    \boldsymbol{\tilde{x}}'=
    \left[\begin{array}{c}
    x'\\
    y'\\
    1
    \end{array}\right]
```

$\boldsymbol{H}$は任意の3×3の行列である。
座標$(x',y')$を求めると、

```math
    x'=\frac{h_{11}x+h_{12}y+h_{13}}{h_{31}x+h_{32}y+h_{33}}\\
    y'=\frac{h_{21}x+h_{22}y+h_{23}}{h_{31}x+h_{32}y+h_{33}}
```

```python
pts1 = np.float32([[0,0], [width,0], [0,height], [width,height]])
pts2 = np.float32([[0,0], [width,50], [0,height], [x1,250]])

#M, _ = cv2.findHomography(pts1, pts2, cv2.RANSAC, 5.0)
M = cv2.getPerspectiveTransform(pts1, pts2)

dst = cv2.warpPerspective(img_rgb, M, (width,height))
plt.imshow(dst)
plt.xticks([]);
plt.yticks([]);
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/1ce0c2e8-8044-8136-e055-08491225ba68.png)


### 合成変換
線形変換の合成変換の性質は、同次座標を利用した射影変換による合成変換に対しても成り立つ。

```math \boldsymbol{\tilde{x}''}\sim(\boldsymbol{H_{2}H_{1}})\boldsymbol{\tilde{x}}
```

## 画像の再標本化と補間
### 画像の再標本化
ディジタル画像に幾何学的変換を施すと、元の画像の画素位置は、一般に変換後は標本化位置からずれる。
変換後の画像を再び縦横等間隔に標本化された位置の値の集まりとして表現するために,**再標本化**が必要となる。
1. 変換後の出力画像におけるある画素位置に対し、適用する幾何学的変換の逆変換を行い、元の入力画像に対する位置を求める。
2. 上で求めた入力画像での位置は、一般に入力画像の画素位置からずれている。そこで、その位置の値を**補間**により周囲の画素値から求める。
3. 上の処理を、出力画像上のすべての画素位置に対して行う。

### ニアレストネイバー
求めたい位置に最も近い画素位置の値をそのまま利用する補間法を、**ニアレストネイバー**と呼ぶ。
求めたい位置を$(x,y)$とすると、その位置の画素値$I(x,y)$は以下のように得られる。

```math
    I(x,y)=f([x+0.5],[y+0.5])
```

ただし、$f(i,j)$は入力画像の位置$(i,j)$の画素値を、記号[ ]はガウス記号で[ ]内の数字を超えない最大の整数値を返すものである。
処理が単純で高速であるが、滑らかなエッジがギザギザになって現れるジャギーが発生しやすい。

### バイリニア補間
求めたい位置$(x,y)$の値$I(x,y)$を、まわりの4点の画素値を用い、以下のように求める。

```math
    I(x,y)=
    \left[\begin{array}{cc}
    [x]+1-x&x-[x]
    \end{array}\right]
    \left[\begin{array}{cc}
    f([x],[y])&f([x],[y]+1)\\
    f([x]+1,[y])&f([x]+1,[y]+1)
    \end{array}\right]
    \left[\begin{array}{c}
    [y]+1-y\\
    y-[y]
    \end{array}\right]\\
    =([x]+1-x)([y]+1-y)f([x],[y])+([x]+1-x)(y-[y])f([x],[y]+1)\\
    +([x]-x)([y]+1-y)f([x]+1,[y])+(x-[x])(y-[y])f([x]+1,[y]+1)
```

これを**バイリニア補間**と呼ぶ。
まわりの4点の画素値の重み付きの平均値を求めることになるため、平滑化の効果が生じる。
そのため、ニアレストネイバーのようなジャギーが目立たなくなるが、エッジがなまってしまう傾向がある。

### バイキュービック補間
**バイキュービック補間**では、求めたい位置$(x,y)$の値$I(x,y)$を、まわりの16点の画素値$f_{11},f_{12},\cdots,f_{44}$を用いて、以下の式により求める。

```math
    I(x,y)=
    \left[\begin{array}{cccc}
    h(x_{1})&h(x_{2})&h(x_{3})&h(x_{4})
    \end{array}\right]
    \left[\begin{array}{cccc}
    f_{11}&f_{12}&f_{13}&f_{14}\\
    f_{21}&f_{22}&f_{23}&f_{24}\\
    f_{31}&f_{32}&f_{33}&f_{34}\\
    f_{41}&f_{42}&f_{43}&f_{44}
    \end{array}\right]
    \left[\begin{array}{c}
    h(y_{1})\\
    h(y_{2})\\
    h(y_{3})\\
    h(y_{4})
    \end{array}\right]\\
```

ただし、

```math
    x_{1}=1+x-[x]\\
    x_{2}=x-[x]\\
    x_{3}=[x]+1-x\\
    x_{4}=[x]+2-x\\
    y_{1}=1+y-[y]\\
    y_{2}=y-[y]\\
    y_{3}=[y]+1-y\\
    y_{4}=[y]+2-y
```

関数$h(t)$は、sinc関数を3次多項式で近似するもので、一般に以下の式が用いられる。

```math
    h(t)=
    \begin{cases}
        |t|^3-2|t|^2+1 & (|t|\leqq1) \\
        -|t|^3+5|t|^2-8|t|+4 & (1<|t|\leqq2) \\
        0 & (2<|t|)
    \end{cases}
```

バイリニア補間に比べて、よりシャープで自然な画像が得られる。

```python
pols = [cv2.INTER_NEAREST, cv2.INTER_LINEAR, cv2.INTER_CUBIC]
pol_names = ['INTER_NEAREST', 'INTER_LINEAR', 'INTER_CUBIC']

fig, ax = plt.subplots(1, 3, figsize=(16, 5))
for i in range(3):
    tmp = cv2.resize(img_rgb, dsize=None, fx=10, fy=10, interpolation = pols[i])
    ax[i].imshow(tmp[500:1000, 500:1000,:])
    ax[i].set_xticks([])
    ax[i].set_yticks([])
    ax[i].set_title(pol_names[i])
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/de9800c1-d714-14a8-095d-aa06e6e8977d.png)


## イメージモザイキング
### イメージモザイキングとその概略処理手順
複数の画像をつなぎ合わせて1つの画像とすることを、**イメージモザイキング**と呼ぶ。
概略の処理手順は、以下のようになる。
1. 特徴点と検出のマッチング
2. 幾何学的変換の推定
3. 画像の幾何学的変換と合成

### 特徴点の検出とマッチング
各入力画像に対し、各々から特徴点を検出する。つぎに、検出された特徴点に対し、画像間で対応を求める。
このように、画像間で対応を求めることを**マッチング**と呼ぶ。

### 幾何学的変換の推定
イメージモザイキングでは射影変換が良く利用される。

```math
    \left[\begin{array}{c}
    x'\\
    y'\\
    1
    \end{array}\right] \sim
    \left[\begin{array}{ccc}
    h_{11}&h_{12}&h_{13}\\
    h_{21}&h_{22}&h_{23}\\
    h_{31}&h_{32}&h_{33}
    \end{array}\right]
    \left[\begin{array}{c}
    x\\
    y\\
    1
    \end{array}\right]
```

は、定数倍を許した方程式であるため$h_{33}=1$とおき、整理することで以下の式を得ることができる。

```math
    xh_{11}+yh_{12}+h_{13}-xx'h_{31}-x'yh_{32}=x'\\
    xh_{21}+yh_{22}+h_{23}-xy'h_{31}-yy'h_{32}=y'
```

未知数の数は8個なので、対応する特徴点が4組以上あれば解を得ることができる。
まず、各対応点から得られる式を以下に示すように並べて行列表現する。

```math
    \boldsymbol{Ah}=\boldsymbol{b}
```

ここで、

```math
    \boldsymbol{A}=
    \left[\begin{array}{cccccccc}
    x_{1}&y_{1}&1&0&0&0&-x_{1}x'_{1}&-x'_{1}y_{1}\\
    0&0&0&x_{1}&y_{1}&1&-x_{1}y'_{1}&-y_{1}y'_{1}\\
    x_{2}&y_{2}&1&0&0&0&-x_{2}x'_{2}&-x'_{2}y_{2}\\
    0&0&0&x_{2}&y_{2}&1&-x_{2}y'_{2}&-y_{2}y'_{2}\\
    &\vdots&&\vdots&&\vdots&&\vdots&
    \end{array}\right],
    \boldsymbol{h}=
    \left[\begin{array}{c}
    h_{11}\\
    h_{12}\\
    h_{13}\\
    h_{21}\\
    h_{22}\\
    h_{23}\\
    h_{31}\\
    h_{32}
    \end{array}\right],
    \boldsymbol{b}=
    \left[\begin{array}{c}
    x'_{1}\\
    y'_{1}\\
    x'_{2}\\
    y'_{2}\\
    \vdots
    \end{array}\right]
```

最小二乗解は、以下の式で得られる。

```math
\boldsymbol{h}=(\boldsymbol{A}^{T}\boldsymbol{A})^{-1}\boldsymbol{A}^{T}\boldsymbol{b}

```

マッチングの結果には、一般には**アウトライヤ**と呼ばれる誤った対応が多数含まれる。
アウトライヤに影響されずに何らかのパラメータを推定する方法として、**RANSAC(RAndom SAmple Conensus)**と呼ばれる方法がよく利用される。
RNSACによるパラメータの推定の手順を、以下に述べる。  
1. パラメータを推定するために必要最小限のデータをランダムに選択する。上の場合、4組の対応点を選択する。
2. 選択した4組の対応点座標を用いて、射影変換のパラメータ$h_{11},h_{12},\dots,h_{32}$を求める。
3. 求めた射影変換パラメータを用いて、他の全ての対応点が正しく変換されているかチェックする。
ここで、正しく変換された対応点は**インライヤ**と呼ばれる。用いたパラメータに対するインライヤの数を記憶する。
4. 1~3の処理を繰り返し、最後に最もインライヤの数が多い射影変換パラメータとそのインライヤを結果として出力する。


画像を2つに分けて片方を変形してつなぎあわせてみます。

```python
img1 = img_rgb[0:400, 0:400]
img2 = img_rgb[50:350, 300:700]

x1 = img2.shape[1]
x2 = img2.shape[0]

pts1 = np.float32([[0,0], [x1,0], [0,x2], [x1,x2]])
pts2 = np.float32([[0,0], [x1,30], [0,x2], [x1,280]])

#M, _ = cv2.findHomography(pts1, pts2, cv2.RANSAC, 5.0)
M = cv2.getPerspectiveTransform(pts1, pts2)

img2 = cv2.warpPerspective(img2, M, (x1, x2))

fig, ax = plt.subplots(1, 2, figsize=(8, 3))

ax[0].imshow(img1);
ax[0].set_xticks([]);
ax[0].set_yticks([]);
ax[1].imshow(img2);
ax[1].set_xticks([]);
ax[1].set_yticks([]);
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/a0979dde-7e57-4a15-f280-bcb6d8a6143a.png)

ここで特徴点のマッチングなどを行いますが、これは次回以降で説明をします。

```python
detector = cv2.AKAZE_create()
matcher = cv2.DescriptorMatcher_create("BruteForce-Hamming")

k1, d1 = detector.detectAndCompute(img1,None)
k2, d2 = detector.detectAndCompute(img2,None)
 
match = cv2.BFMatcher()
matches = match.knnMatch(d2, d1, k=2)

good = []
for m, n in matches:
    if m.distance < 0.8* n.distance:
        good.append(m)

MIN_MATCH_COUNT = 10
if len(good) > MIN_MATCH_COUNT:
    src_pts = np.float32([ k2[m.queryIdx].pt for m in good ])
    dst_pts = np.float32([ k1[m.trainIdx].pt for m in good ])
    H = cv2.findHomography(src_pts, dst_pts, cv2.RANSAC, 5.0)[0]
else:
    print('Not enought matches are found - {}/{}'.format(len(good), MIN_MATCH_COUNT))
    exit(1)

img2_warped = cv2.warpPerspective(img2, H, (img1.shape[1] + img2.shape[1], img1.shape[0]))

img_stitched = img2_warped.copy()
img_stitched[:img1.shape[0], :img1.shape[1]] = img1

def trim(frame):
    if np.sum(frame[0]) == 0:
        return trim(frame[1:])
    if np.sum(frame[-1]) == 0:
        return trim(frame[:-2])
    if np.sum(frame[:,0]) == 0:
        return trim(frame[:, 1:])
    if np.sum(frame[:,-1]) == 0:
        return trim(frame[:, :-2])
    return frame
img_stitched_trimmed = trim(img_stitched)

img_key = cv2.drawKeypoints(img2, k2, None)
draw_params = dict(matchColor=(0,255,0),
                   singlePointColor=None,
                   flags=2)
img_match = cv2.drawMatches(img2, k2, img1, k1, good, None, **draw_params)
```

```python
plt.imshow(img_match)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/11ff7acf-5121-28c5-f231-3f3b02ad75e2.png)

```python
plt.imshow(img_stitched_trimmed)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/24a17bdf-c253-9e23-7b3c-26252d456b90.png)


### 画像の幾何学的変換と合成
画像をつなぎ合わせたとき、画像のつなぎ目が目立ってしまうことがある。
画像のつなぎ目の位置を適切に選ぶか、重なった部分で両方の画素値を混ぜあわせる方法がある。
後者の例として、両方の重なった部分でアルファブレンディングの考え方を利用し、画素値を重み付き平均する方法がある。
例えば、画像1の画素値を$I_{1}$、画像2の画素値を$I_{2}$、画像1の画像端からの距離を$d_{1}$、画像2の画像端からの距離を$d_{2}$とするとき、
以下の式により、つなぎ目を目立たなくする。

```math
    \frac{d_{1}I_{1}+d_{2}I_{2}}{d_{1}+d_{2}}
```

これで幾何学変換についての勉強は終わります。

### 次回
2値画像処理

### 参考

https://www.amazon.co.jp/%E3%83%87%E3%82%A3%E3%82%B8%E3%82%BF%E3%83%AB%E7%94%BB%E5%83%8F%E5%87%A6%E7%90%86-%E6%94%B9%E8%A8%82%E7%AC%AC%E4%BA%8C%E7%89%88-%E3%83%87%E3%82%A3%E3%82%B8%E3%82%BF%E3%83%AB%E7%94%BB%E5%83%8F%E5%87%A6%E7%90%86%E7%B7%A8%E9%9B%86%E5%A7%94%E5%93%A1%E4%BC%9A/dp/490347464X