以前、台形補正についての記事を書きました。

今回は、台形だけではなく、任意の四角形を正方形（または別の任意の四角形）に変換しよう、という試み。

射影変換の原理を勉強して、Pythonで実装してみました。

PythonライブラリのSympyとScikit-imageを使う方法です。

1. 射影変換の原理

クリックすると展開されます

射影変換（ホモグラフィ）の原理にはついて、たくさんの解説記事が見つかります。

shogo82148.github.io

yaju3d.hatenablog.jp

この記事では、私が自分なりに理解したことについて書いていきます。

・ステップ1

ベクトル $(x,y)$ をベクトル $(u,v)$ にうつす線形写像 $H$ を求める問題を考えます。

f:id:plant-raspberrypi3:20181118185617p:plain:w400

これは

${\displaystyle \mathbf{a} = \begin{pmatrix} x \\ y \\ 1 \end{pmatrix} \\ }$

${\displaystyle \mathbf{b} = w \begin{pmatrix} u \\ v \\ 1 \end{pmatrix} = \begin{pmatrix} uw \\ vw \\ w \end{pmatrix} \\ }$

${\displaystyle H = \begin{pmatrix} h_{11} & h_{12} & h_{13} \\ h_{21} & h_{22} & h_{23} \\ h_{31} & h_{32} & 1 \end{pmatrix} }$

とおくと

$\mathbf{b} = H \mathbf{a}$

が成り立つ $H$ を求める問題と置き換えることができます。

・ステップ2

次に、四隅の座標を基準にして変形する場合を考えます。変換前の座標を

$(x_{n}, y_{n})$

変換後の座標を

$(u_{n}, v_{n})$

とします。ただし、 $n = 1, 2, 3, 4$ です。

具体的には、以下の図のような変換を考えます。

f:id:plant-raspberrypi3:20181118190217p:plain:w400

ここで、四隅の座標が与えられている場合、各 $n$ について以下の式が成り立ちます。

$u_{n} w_{n} = h_{11} x_{n} + h_{12} y_{n} + h_{13}$
$v_{n} w_{n} = h_{21} x_{n} + h_{22} y_{n} + h_{23}$

これは、 $\mathbf{b} = H \mathbf{a}$ を計算することで得られます。

また、同様に

$w_{n} = h_{31} x_{n} + h_{32} y_{n} + 1$

も得られます。これを上の2つの式に代入して左辺を右辺に移動すると

$h_{11} x_{n} + h_{12} y_{n} + h_{13} - u_{n} (h_{31} x_{n} + h_{32} y_{n} + 1) = 0$
$h_{21} x_{n} + h_{22} y_{n} + h_{23} - v_{n} (h_{31} x_{n} + h_{32} y_{n} + 1) = 0$

となるので、これを $n = 1, 2, 3, 4$ について解くことで $h_{11}$ 〜 $h_{32}$ が得られ、 $H$ が求まります。

求めた $H$ を使い、

$\mathbf{b} = H\mathbf{a}$

の関係から画像の変換を行います。

2. Sympyで連立方程式を解きHを求める

クリックすると展開されます

前述の通り、 $H$ を求めるためには8元連立方程式を解く必要があります。

しかし、手計算で解くのはとても面倒なので、Sympyで解いていきます。

Sympyの使い方についてはこの辺りを参照。

3.2. Sympy : Python での代数計算 — Scipy lecture notes

環境は、MacOS Sierra, Python3.6です。まずは必要なパッケージをインポートします。

from sympy import solve,symbols
import numpy as np

$h_{11}$ 〜 $h_{32}$ を変数として設定します。

str_list = [f"h{i}" for i in [11,12,13,21,22,23,31,32]]
str_symbol = " ".join(str_list)

h11, h12, h13, h21, h22, h23, h31, h32 = symbols(str_symbol)

次に、 $(x_{n}, y_{n})$ のセットをnumpyアレイxyに、 $(u_{n}, v_{n})$ のセットをnumpyアレイuvに代入します。

$H$ を解くための連立方程式を以下のように記述します。

各 $n$ ごとに方程式を求め、リストexpに加えます。

exp = []
for i in range(4):
    exp_x = h11*xy[i,0]+h12*xy[i,1]+h13-(h31*xy[i,0]+h32*xy[i,1]+1)*uv[i,0]
    exp_y = h21*xy[i,0]+h22*xy[i,1]+h23-(h31*xy[i,0]+h32*xy[i,1]+1)*uv[i,1]
    exp.append(exp_x)
    exp.append(exp_y)

全ての方程式の入ったリストexpを引数としてSympyのsolve関数に渡し、連立方程式を解き、解をHに代入します。

result = solve(exp,dict=True)
r = result[0]

H = np.array([[r[h11],r[h12],r[h13]],
                           [r[h21],r[h22],r[h23]],
                           [r[h31],r[h32],1]],dtype=np.float64)

仮に以下のような図の変換の場合

f:id:plant-raspberrypi3:20181118194300p:plain

方程式は

1: 10.0*h11 + 20.0*h12 + h13
2: 10.0*h21 + 20.0*h22 + h23
3: 20.0*h11 + 70.0*h12 + h13
4: 20.0*h21 + 70.0*h22 + h23 - 2000.0*h31 - 7000.0*h32 - 100.0
5: 90.0*h11 + 90.0*h12 + h13 - 9000.0*h31 - 9000.0*h32 - 100.0
6: 90.0*h21 + 90.0*h22 + h23 - 9000.0*h31 - 9000.0*h32 - 100.0
7: 90.0*h11 + 10.0*h12 + h13 - 9000.0*h31 - 1000.0*h32 - 100.0
8: 90.0*h21 + 10.0*h22 + h23

となり、Sympyによって $H$ を求めると

array([[ 2.06775593e+00, -4.13551186e-01, -1.24065356e+01],
       [ 2.13237330e-01,  1.70589864e+00, -3.62503462e+01],
       [ 8.18794424e-03, -4.13551186e-03,  1.00000000e+00]])

となることがわかりました。

この $H$ が正しいかは、 $(x_{n}, y_{n})$ の値を用いて $H\mathbf{a}$ を計算して確かめることができます。

for i in range(4):
    xy1 = np.append(xy[i],1)
    uv1 = np.dot(H,xy1)
    uv1 = np.round(uv1/uv1[-1]).astype(int)
    uv1 = uv1[:-1]
    print(f"{xy[i]} -> {uv1}")

#出力
[20 10] -> [0 0]
[70 20] -> [100   0]
[90 90] -> [100 100]
[10 90] -> [  0 100]

3. 求めたHとScikit-imageを使って図形を変換する

クリックすると展開されます

$H$ を図形に適用して、実際に変換していきます。

変換がわかりやすいように多角形を描いてみました。

polygon = np.array([[20,10],[50,40],[70,20],[90,90],[50,60],[10,90]])

img = np.zeros((100,100),dtype=np.float64)
rr,cc = draw.polygon(polygon.T[0],polygon.T[1])
img[rr,cc] = 1

plt.imshow(img)
plt.show()

f:id:plant-raspberrypi3:20181118203324p:plain:w400

これを変換していきます。

前の項目で述べた、 $H$ を求める過程を関数にまとめます。

from sympy import solve,symbols
import numpy as np
from skimage import draw,transform
import matplotlib.pyplot as plt
from pprint import pprint

def solve_H(xy,uv):

    print("XY to UV")
    for i in range(4):
        print(f"{xy[i]} -> {uv[i]}")

    str_list = [f"h{i}" for i in [11,12,13,21,22,23,31,32]]
    str_symbol = " ".join(str_list)
    
    h11, h12, h13, h21, h22, h23, h31, h32 = symbols(str_symbol)

    exp = []

    for i in range(4):
        exp_x = h11*xy[i,0]+h12*xy[i,1]+h13-(h31*xy[i,0]+h32*xy[i,1]+1)*uv[i,0]
        exp_y = h21*xy[i,0]+h22*xy[i,1]+h23-(h31*xy[i,0]+h32*xy[i,1]+1)*uv[i,1]
        exp.append(exp_x)
        exp.append(exp_y)

    print("")
    print("Equations")
    for i,e in enumerate(exp):
        print(f"{i+1}: {e}")

    result = solve(exp,dict=True)
    r = result[0]

    H = np.array([[r[h11],r[h12],r[h13]],
                           [r[h21],r[h22],r[h23]],
                           [r[h31],r[h32],1]],dtype=np.float64)
    print("")
    print("Result_array")
    pprint(H)
    
    return H

この関数を使って解きます。四隅の変換前と変換後の座標を設定します。

xy = np.array([[20,10],[70,20],[90,90],[10,90]])
uv = np.array([[0,0],[100,0],[100,100],[0,100]])

早速、 $H$ を解いていきます。

が、そのままでやると何故かわからないですがうまくいかないことがわかったので、xとyを入れ替えたりゴソゴソやっています。numpyのxyと画像にした時のxyの方向が異なるせいかなぁ、と思っていますが、もし詳しい方がいらっしゃったらコメントいただけるとありがたいです。。。

xy_t = xy[:,::-1]
uv_t = uv[:,::-1]

H = solve_H(xy_t,uv_t)

#出力
XY to UV
[10 20] -> [0 0]
[20 70] -> [  0 100]
[90 90] -> [100 100]
[90 10] -> [100   0]

Equations
1: 10.0*h11 + 20.0*h12 + h13
2: 10.0*h21 + 20.0*h22 + h23
3: 20.0*h11 + 70.0*h12 + h13
4: 20.0*h21 + 70.0*h22 + h23 - 2000.0*h31 - 7000.0*h32 - 100.0
5: 90.0*h11 + 90.0*h12 + h13 - 9000.0*h31 - 9000.0*h32 - 100.0
6: 90.0*h21 + 90.0*h22 + h23 - 9000.0*h31 - 9000.0*h32 - 100.0
7: 90.0*h11 + 10.0*h12 + h13 - 9000.0*h31 - 1000.0*h32 - 100.0
8: 90.0*h21 + 10.0*h22 + h23

Result_array
array([[ 2.06775593e+00, -4.13551186e-01, -1.24065356e+01],
       [ 2.13237330e-01,  1.70589864e+00, -3.62503462e+01],
       [ 8.18794424e-03, -4.13551186e-03,  1.00000000e+00]])

最後に求めた $H$ を、skimageのtransformモジュールの機能を使って、画像に適用します。
使う機能はProjectiveTransformとwarpです。

ProjectiveTransformはinverseにすると想定通りの変換となります。

変換によって伸展した場所がモヤっとしてしまうので、最後に四捨五入によってバイナリ画像に修正してから表示します。

tform = transform.ProjectiveTransform(matrix=H).inverse
warped = transform.warp(img,tform,order=3,mode="edge")

plt.imshow(img)
plt.show()

plt.imshow(np.round(warped).astype(int))
plt.show()

変換前

f:id:plant-raspberrypi3:20181118200712p:plain:w400