行列によるアフィン変換（拡大縮小・回転・剪断・移動）　～Python画像処理の再発明家～

画像処理ライブラリに頼らず、行列演算だけでアフィン変換をするお話。Pythonistaでも可能

**[基礎編はこちらから](http://qiita.com/secang0/items/1229212a37d8c9922901)**

#「再発明家」とは
Open CVとかPillowに頼らず、numpyとmatplotlibを使って、様々な画像処理を実際に書いてみる。iOSアプリのPythonistaでも使える組み合わせだ。

```py3
import numpy as np
import matplotlib.pyplot as plt
```

また、画像表示には以下の関数が便利である。（詳しくは[基礎編](http://qiita.com/secang0/items/1229212a37d8c9922901)）

```py3
def img_show(img : np.ndarray, cmap = 'gray', vmin = 0, vmax = 255, interpolation = 'none') -> None:
    '''np.arrayを引数とし、画像を表示する。'''
    
    #dtypeをuint8にする
    img = np.clip(img,vmin,vmax).astype(np.uint8)
    
    #画像を表示
    plt.imshow(img, cmap = cmap, vmin = vmin, vmax = vmax, interpolation = interpolation)
    plt.show()
    plt.close()
```

#アフィン変換と行列
画像をゆがませる方法はさまざまある。線形変換（拡大縮小・回転・剪断）と平行移動を合わせた変換をアフィン変換という。
[ここ](http://d.hatena.ne.jp/Zellij/20120523/p1)が分かりやすかった。ググっても分かりやすい説明はいくらでも転がっている。

このアフィン変換は行列の掛け算で表現できる。
あるアフィン変換$A$が点$(x_0,y_0)$から$(x_1,y_1)$に移したとき、

```math
\left(
\begin{matrix}
x_1\\
y_1\\
1
\end{matrix}
\right)
=A
\left(
\begin{matrix}
x_0\\
y_0\\
1
\end{matrix}
\right)
```

と書くことができる。この$A$の係数を見てみると

```math
\left(
\begin{matrix}
a &b& t_x\\
c &d&y_y\\
0&0&1
\end{matrix}
\right)
```

という形をしている。このうち、$a,b,c,d$は変形を、$t_x,t_y$は平行移動を担当する。
一応計算しておくと、

```math
\left(
\begin{matrix}
x_1\\
y_1\\
1
\end{matrix}
\right)
=\left(
\begin{matrix}
a &b& t_x\\
c &d&t_y\\
0&0&1
\end{matrix}
\right)
\left(
\begin{matrix}
x_0\\y_0\\1
\end{matrix}
\right)
=
\left(
\begin{matrix}
ax_0 +by_0+ t_x\\
cx_0 +dy_0+t_y\\
0+0+1
\end{matrix}
\right)
```

さて、画像変換でこれを応用してみよう。
['tiger.jpeg'](http://imgur.com/AfaVcBI)を切り出して使った。

```py3
img = plt.imread('tiger.jpeg')[1390:1440,375:425]
img_show(img)
```

![tiger_trimmed_E.png](https://qiita-image-store.s3.amazonaws.com/0/172464/23f0d49f-c26e-35cb-d257-a19130bdc983.png)


手順は
1. 変換後の各ピクセルの座標を計算
2. 変換後の各ピクセルが参照すべき、変換前の座標を計算
3. 変換前の座標に応じて、各ピクセルの値を決定


#アフィン変換後の各ピクセルの座標を計算

まずは、各ピクセルについて、その座標を収めた二次元配列を作ろう。上のベクトルを真似て、ついでに1を末尾につけておく。

```py3
#高さ3、幅4の画像を作る
height, width = 3,4

#mgridでx座標の行列、y座標の行列を作成
x, y = np.mgrid[:x_len,:y_len]

#dstackでx座標、y座標、1を組み合わせる
xy_after = np.dstack((x,y,np.ones((x_len, y_len))))
xy_after

#array([
#[[ 0.,  0.,  1.], [ 0.,  1.,  1.], [ 0.,  2.,  1.], [ 0.,  3.,  1.]],
#[[ 1.,  0.,  1.], [ 1.,  1.,  1.], [ 1.,  2.,  1.], [ 1.,  3.,  1.]],
#[[ 2.,  0.,  1.], [ 2.,  1.,  1.], [ 2.,  2.,  1.], [ 2.,  3.,  1.]]])
```


#変換後の各ピクセルが参照すべき、変換前の座標を計算
画像処理では、アフィン行列を直接用いるのではなく、その逆行列を用いる（正則であると考える）
その理由は「参照すべき座標を各ピクセルについて決定するためである」

```py3
#縦横に2倍拡大するアフィン変換
affin = np.matrix('2,0,0;0,2,0;0,0,1')

#逆行列
inv_affin = np.linalg.inv(affin)

#行列の掛け算をアインシュタイン和で計算
ref_xy = np.einsum('ijk,lk->ijl',xy_after,inv_affin)[...,:2]
ref_xy

#array([
#[[ 0. ,  0. ], [ 0.5,  0. ], [ 1. ,  0. ]],
#[[ 0. ,  0.5], [ 0.5,  0.5], [ 1. ,  0.5]],
#[[ 0. ,  1. ], [ 0.5,  1. ], [ 1. ,  1. ]],
#[[ 0. ,  1.5], [ 0.5,  1.5], [ 1. ,  1.5]]])
```

このように、例えば変換後の$(1,1)$が変換前の$(0.5,0.5)$と一致する、と知るために逆行列を用いたのである。

#変換前の座標に応じて、各ピクセルの値を決定
上の行列'ref_xy'を見ると`[2.,2.]`は`[1.,1.]`のピクセルの値と一致させればよいことが分かる。しかし、`[1.,2.]`などは`[0.5,1.]`という存在しないピクセルを参照しなければいけない。この存在しない座標をどうするか。

以下その方法を2つ紹介したい。
[ここを見るとわかりやすい](http://ipr20.cs.ehime-u.ac.jp/column/gazo_syori/chapter3.html)

##最近傍法
簡単に言うと四捨五入を行う手法である。四捨五入には0.5を足してint型に変換すればよい。

以下のコードは実際に画像を拡大している。

```py3
#参照する座標が四捨五入で計算されるため、100,450にするとインデックスエラーになる
height, width = 99, 149

x,y = np.mgrid[:height,:width]
xy_after = np.dstack((x,y,np.ones((height, width))))

#アフィン変換の行列を用意
#縦に2倍、横に3倍
affin = np.matrix('2,0,0;0,3,0;0,0,1')
inv_affin = np.linalg.inv(affin)

#参照する座標を計算
ref_xy = np.einsum('ijk,lk->ijl',xy_after,inv_affin)[...,:2]
ref_nearmost_xy = (ref_xy + 0.5).astype(int)
img_nearmost = img[ref_nearmost_xy[...,0],ref_nearmost_xy[...,1]]

img_show(img_nearmost)
```
![nearmost.png](https://qiita-image-store.s3.amazonaws.com/0/172464/322ee126-fd61-7f31-c065-22c8af83dc3f.png)



#線形補間法
[先ほどのリンク](http://ipr20.cs.ehime-u.ac.jp/column/gazo_syori/chapter3.html)の写真がやはりわかりやすい。
![](http://ipr20.cs.ehime-u.ac.jp/column/gazo_syori/image/zahyou.gif)

この方法では、近い4つのピクセルを、そのピクセルの近さで重みづけしている。

まず、近いピクセルを計算する。

```py3
#左上をintで計算した後、それを移動させて計算
linear_xy = {}
linear_xy['upleft'] = ref_xy.astype(int)
linear_xy['downleft'] = linear_xy['upleft'] + [1,0]
linear_xy['upright']= linear_xy['upleft'] + [0,1]
linear_xy['downright'] = linear_xy['upleft'] + [1,1]
```
次に、左上のピクセルとの差を計算して、重みづけを計算する。

```py3
#左上の点との差を計算
upleft_diff = ref_xy - linear_xy['upleft']

#(1-xの差)と(1-yの差)の積を計算
linear_weight = {}
linear_weight['upleft'] = (1-upleft_diff[...,0])*(1-upleft_diff[...,1])
linear_weight['downleft'] = upleft_diff[...,0]*(1-upleft_diff[...,1])
linear_weight['upright'] = (1-upleft_diff[...,0])*upleft_diff[...,1]
linear_weight['downright'] = upleft_diff[...,0]*upleft_diff[...,1]
```

あとは、これを掛け合わせて画素値を計算するだけである。

```py3
#height, width = 98, 147
#affin = np.matrix('2,0,0;0,3,0;0,0,1')
#としている

linear_with_weight = {}
for direction in liner_xy.keys():
    xy = linear_xy[direction]
    weight = linear_weight[direction]
    linear_with_weight[direction] = np.einsum('ij,ijk->ijk',weight,img[xy[...,0],xy[...,1]])
img_linear = sum(linear_with_weight.values())
img_show(img_linear)
```

![linear.png](https://qiita-image-store.s3.amazonaws.com/0/172464/617a715a-0467-7592-2f8f-8c153124e1bb.png)

微妙な違いがあり、こっちのほうがなめらかである。

#背景を設定
画像の変形の方法や変形後の形によってはインデックスエラーが出る。
その理由は、存在しない画素を参照しているからである。
とりあえず座標が0を下回ったり、最大値を上回ったもの-1で置き換える関数を定義する。

```py3
def clip_xy(ref_xy, img_shape):
    #x座標について置換
    ref_x = np.where((0<=ref_xy[...,0])&(ref_xy[...,0]<img_shape[0]),ref_xy[...,0],-1)
    #y座標について置換
    ref_y = np.where((0<=ref_xy[...,1])&(ref_xy[...,1]<img_shape[1]),ref_xy[...,1],-1)
    
    #結合して返す
    return np.dstack([ref_x,ref_y])
```

すると実は、-1に置き換えることで、今まで外れたピクセルを参照していたピクセルは、すべて最終行、最終列を参照するようになった。（-1でなく、`img_shape[0]`などでも問題ない）
あとは背景色を持った最終行、最終列を作るだけだ。


```py3
#背景色を設定
bg_color = [0,0,0]

#背景色で埋められた一回り大きい画像を作成
img_bg = np.empty(np.add(img.shape,(1,1,0)))
img_bg[:,:] = bg_color

#画像を貼り付ける
img_bg[:-1,:-1] = img


#高さ150・幅500の変換後画像を作成
height, width = 150, 500
x,y = np.mgrid[:height,:width]
xy_after = np.dstack((x,y,np.ones((height, width))))
#アフィン変換の行列を用意
#縦に2倍、横に3倍
affin = np.matrix('2,0,0;0,3,0;0,0,1')
inv_affin = np.linalg.inv(affin)

#最近傍法で画像を変換
ref_xy = np.einsum('ijk,lk->ijl',xy_after,inv_affin)[...,:2]
ref_nearmost_xy = (ref_xy + 0.5).astype(int)
ref_nearmost_xy = clip_xy(ref_nearmost_xy)
#clip_xyで参照するピクセルが最終行、最終列を参照するように変更
img_nearmost_bg = img_bg[ref_nearmost_xy[...,0],ref_nearmost_xy[...,1]]
img_show(img_nearmost_bg)
```

![nearmost_bg.png](https://qiita-image-store.s3.amazonaws.com/0/172464/d78dd1e7-b844-ae2d-3ca2-22191f4b9636.png)

このように黒色の背景が追加されている。

あとはアフィン変換の行列を変更して、自由に遊んでほしい。

`affin = np.matrix('2,0.5,15;1,-3,200;0,0,1')`
![free.png](https://qiita-image-store.s3.amazonaws.com/0/172464/60900444-04d3-9c09-0c12-dfbf8216e8f4.png)