つくるって楽しい

主にpythonとか。画像処理とか。

移動座標の対応点からAffine変換の係数を求める

3組の対応点からAffine変換を求める方法について。opencvだとcv2.getAffineTransform()に対応します。なぜ3点なのでしょう。

理論

Affine変換により (x_{i}, y_{i}) (x_{i}^{'}, y_{i}^{'})に変換される時次のように書けます。

{
\begin{pmatrix}
x_{i}^{'} \\
y_{i}^{'} \\
\end{pmatrix}

=

\begin{pmatrix}
a_{11} & a_{12} & a_{13} \\
a_{21} & a_{22} & a_{23} \\
\end{pmatrix}

\begin{pmatrix}
x_{i} \\
y_{i} \\
1 \\
\end{pmatrix}
\tag{1}
}

(1)を書き換えます。

{
\displaystyle
\begin{pmatrix}
x_{i} && y_{i} && 1 && 0 && 0 && 0 \\
0 && 0 && 0 && x_{i} && y_{i} && 1 \\
\end{pmatrix}

\begin{pmatrix}
a_{11} \\
a_{12} \\
a_{13} \\
a_{21} \\
a_{22} \\
a_{23} \\
\end{pmatrix}

=

\begin{pmatrix}
x_{i}^{'} \\
y_{i}^{'} \\
\end{pmatrix}
\tag{2}
}

(2)は1対応点のみの場合ですが、N個の対応点がある場合は次のように拡張できます。

{
\displaystyle
\begin{pmatrix}
x_{1} && y_{1} && 1 && 0 && 0 && 0 \\
0 && 0 && 0 && x_{1} && y_{1} && 1 \\
\vdots &&  &&  &&  &&  && \vdots \\
x_{N} && y_{N} && 1 && 0 && 0 && 0 \\
0 && 0 && 0 && x_{N} && y_{N} && 1 \\
\end{pmatrix}

\begin{pmatrix}
a_{11} \\
a_{12} \\
a_{13} \\
a_{21} \\
a_{22} \\
a_{23} \\
\end{pmatrix}

=

\begin{pmatrix}
x_{1}^{'} \\
y_{1}^{'} \\
\vdots \\
x_{N}^{'} \\
y_{N}^{'} \\
\end{pmatrix}
\tag{3}
}

(3)を以下のように表記します。
{
X \boldsymbol{a} = X^{'}
\tag{4}
}

これは最小二乗法の行列表現[1]ですね。Xが正則ならば\boldsymbol{a}
{
\boldsymbol{a} = X^{-1} X^{'}
\tag{5}
}
と求めることが出来ます。\boldsymbol{a}の未知数は6つなので、対応点が3つあれば(5)の方法で求めることが出来ます。

また、対応点が4つ以上の場合はXが正則でないので逆行列が求められませんが、以下の式で解くことが出来ます。
{
\boldsymbol{a} = (X^{T} X)^{-1} X^{T} X^{'}
\tag{6}
}

(6)は(X^{T} X)^{-1}が正則な時の唯一つの解になるそうです。[1]

試してみた

時計回りに60度、x,y方向にそれぞれ5,3ずつ移動する変換Mを用意します。((1)の形式の変換行列です。)

theta = 60
tx, ty = 5, 3
M = np.array([[np.cos(np.deg2rad(theta)), -np.sin(np.deg2rad(theta)), tx],
                 [np.sin(np.deg2rad(theta)), np.cos(np.deg2rad(theta)), ty]], float)

Mを出力。

M: 
[[ 0.5       -0.8660254  5.       ]
 [ 0.8660254  0.5        3.       ]]

3点の場合

座標変換前の点の集合。

X = np.array([[15, 5, 1],
              [25, 10, 1],
              [20, 15, 1]])

変換。

X_tfm = M.dot(X.T)

前後をプロット。赤:移動前、青:移動後。
f:id:mikekochang:20190227164507j:plain

(5)から変換行列を計算。Mと一致しています。

M_slv = np.linalg.inv(XX).dot(XX_tfm)
M_slv: 
[ 0.5       -0.8660254  5.         0.8660254  0.5        3.       ]

4点の場合

座標変換前の点の集合。変換。

X = np.array([[15, 5, 1],
              [25, 10, 1],
              [20, 15, 1],
              [18, 18, 1]])
X_tfm = M.dot(X.T)

(6)から変換行列を計算。Mと一致しています。

M_slv = np.linalg.inv(XX.T.dot(XX)).dot(XX.T).dot(XX_tfm)
M_slv: 
[ 0.5       -0.8660254  5.         0.8660254  0.5        3.       ]

参考

[1]最小二乗法の行列表現(単回帰,多変数,多項式),https://mathtrain.jp/leastsquarematrix