画像処理ライブラリに頼らず、行列演算だけでアフィン変換をするお話。Pythonistaでも可能
#「再発明家」とは
Open CVとかPillowに頼らず、numpyとmatplotlibを使って、様々な画像処理を実際に書いてみる。iOSアプリのPythonistaでも使える組み合わせだ。
import numpy as np
import matplotlib.pyplot as plt
また、画像表示には以下の関数が便利である。（詳しくは基礎編）
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()
#アフィン変換と行列
画像をゆがませる方法はさまざまある。線形変換（拡大縮小・回転・剪断）と平行移動を合わせた変換をアフィン変換という。
ここが分かりやすかった。ググっても分かりやすい説明はいくらでも転がっている。
このアフィン変換は行列の掛け算で表現できる。
あるアフィン変換$A$が点$(x_0,y_0)$から$(x_1,y_1)$に移したとき、
\left(
\begin{matrix}
x_1\\
y_1\\
1
\end{matrix}
\right)
=A
\left(
\begin{matrix}
x_0\\
y_0\\
1
\end{matrix}
\right)
と書くことができる。この$A$の係数を見てみると
\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$は平行移動を担当する。
一応計算しておくと、
\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'を切り出して使った。
img = plt.imread('tiger.jpeg')[1390:1440,375:425]
img_show(img)
手順は
- 変換後の各ピクセルの座標を計算
- 変換後の各ピクセルが参照すべき、変換前の座標を計算
- 変換前の座標に応じて、各ピクセルの値を決定
#アフィン変換後の各ピクセルの座標を計算
まずは、各ピクセルについて、その座標を収めた二次元配列を作ろう。上のベクトルを真似て、ついでに1を末尾につけておく。
#高さ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.]]])
#変換後の各ピクセルが参照すべき、変換前の座標を計算
画像処理では、アフィン行列を直接用いるのではなく、その逆行列を用いる（正則であると考える）
その理由は「参照すべき座標を各ピクセルについて決定するためである」
#縦横に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つ紹介したい。
ここを見るとわかりやすい
##最近傍法
簡単に言うと四捨五入を行う手法である。四捨五入には0.5を足してint型に変換すればよい。
以下のコードは実際に画像を拡大している。
#参照する座標が四捨五入で計算されるため、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)
#線形補間法
先ほどのリンクの写真がやはりわかりやすい。
この方法では、近い4つのピクセルを、そのピクセルの近さで重みづけしている。
まず、近いピクセルを計算する。
#左上を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]
次に、左上のピクセルとの差を計算して、重みづけを計算する。
#左上の点との差を計算
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]
あとは、これを掛け合わせて画素値を計算するだけである。
#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)
微妙な違いがあり、こっちのほうがなめらかである。
#背景を設定
画像の変形の方法や変形後の形によってはインデックスエラーが出る。
その理由は、存在しない画素を参照しているからである。
とりあえず座標が0を下回ったり、最大値を上回ったもの-1で置き換える関数を定義する。
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]などでも問題ない）
あとは背景色を持った最終行、最終列を作るだけだ。
#背景色を設定
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)
このように黒色の背景が追加されている。
あとはアフィン変換の行列を変更して、自由に遊んでほしい。
affin = np.matrix('2,0.5,15;1,-3,200;0,0,1')
