画像処理ライブラリに頼らず、行列演算だけでアフィン変換をするお話。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. 変換後の各ピクセルの座標を計算
2. 変換後の各ピクセルが参照すべき、変換前の座標を計算
3. 変換前の座標に応じて、各ピクセルの値を決定

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

まずは、各ピクセルについて、その座標を収めた二次元配列を作ろう。上のベクトルを真似て、ついでに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つ紹介したい。
ここを見るとわかりやすい

線形補間法

先ほどのリンクの写真がやはりわかりやすい。

この方法では、近い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')

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