読者です 読者をやめる 読者になる 読者になる

デジタル・デザイン・ラボラトリーな日々

アラフィフプログラマーが数学と物理を基礎からやり直す

【コンピュータビジョン】3章「画像間の写像」 ホモグラフィ(点の調整)

前記事で、「homography.pyでは、更に点群を平均値が0、標準偏差が1になるように正規化していますが、今回はそこまで説明はしません。」と逃げたのだが、アフィン変換でも同様のことをしていたので、頑張って調べてみます。

参考:標準化した値は、平均値がゼロ、分散(と標準偏差)が1になる

H_from_pointsの呼び出し処理

fp = array([[0,120,120,0],[0,0,160,160],[1,1,1,1]])
tp = array([[20,120,120,0],[20,0,160,140],[1,1,1,1]])
H = homography.H_from_points(fp,tp)

H_from_pointsの中身で開始点の部分

# 点を調整する(数値計算上重要)
# 開始点
m = mean(fp[:2], axis=1) #平均
maxstd = max(std(fp[:2], axis=1)) + 1e-9 #標準偏差
C1 = diag([1/maxstd, 1/maxstd, 1]) #対角項
C1[0][2] = -m[0]/maxstd
C1[1][2] = -m[1]/maxstd
fp = dot(C1,fp) #内積

ここから解析していきます。

m = mean(fp[:2], axis=1) #平均

meanは平均で、axis=1:横方向となります。
mには、[0,120,120,0]の平均 60と、[0,0,160,160]の平均 80の2値が返ります。

maxstd = max(std(fp[:2], axis=1)) + 1e-9 #標準偏差

stdは標準偏差で、axis=1:横方向となります。
[0,120,120,0]の標準偏差 60と、と[0,0,160,160]の標準偏差 80の2値が求まり、
maxstdには、60と80の最大値として80に1e-9を加算した80.000000001が返ります。

C1 = diag([1/maxstd, 1/maxstd, 1]) #対角項
C1[0][2] = -m[0]/maxstd
C1[1][2] = -m[1]/maxstd

diagは対角項で、3x3の正方行列を作成します。
1/maxstdの値は0.0125、-m[0]/maxstdの値は-0.75、-m[1]/maxstdの値は-1となります。
\begin{pmatrix}1/maxstd&0&-m[0]/maxstd\\0&1/maxstd&-m[1]/maxstd\\0&0&1\end{pmatrix}=\begin{pmatrix}0.0125&0&-0.75\\0&0.0125&-1\\0&0&1\end{pmatrix}
これは何をしているのか?
アフィン変換に拡大縮小(sx,sy)と平行移動(tx,ty)があります。
\begin{pmatrix}sx&0&tx\\0&sy&ty\\0&0&1\end{pmatrix}
これに当てはめると「1/標準偏差の最大値」で拡大縮小させた後に「平均値/標準偏差の最大値」で平行移動させていると思われます。
この部分が本の記述にあった「平均値が0に、標準偏差が1になるように正規化」にあたるのではないでしょうか。

fp = dot(C1,fp) #内積

dotは内積で下記が計算結果となります。
内積するのは開始点(fp)にアフィン変換行列(C1)を掛けた行列を求めるためです。
参考:文系のための「内積」(2)
\begin{pmatrix}-0.75&0.75&0.75&-0.75\\-1&-1&1&1\\1&1&1&1\end{pmatrix}
この値が何を示しているのを下図で説明します。
この図は、縦幅(120)と横幅(160)から最大値である横幅(160)で正方形にして、赤線が交わったところが重心となります。ちなみに縦幅(120) / 横幅(160)=0.75となります。
f:id:Yaju3D:20140330214806j:plain
開始点の4点座標は左上から反時計周りで(0, 0)、(0,120)、(160,120)、(160, 0)となっています。
内積の計算結果で1行目がY座標、2行目がX座標として開始点と位置を合わせてみます。
重心を原点0として対角線の長さで4点座標とすると下記のようになります。

(0,  0) (160,  0) → (-1,-0.75) (1,-0.75)
(0,120) (160,120) → (-1, 0.75) (1, 0.75)


次はH_from_pointsの中身で対応点の部分

# 対応点
m = mean(tp[:2], axis=1) #平均
maxstd = max(std(fp[:2], axis=1)) + 1e-9 #標準偏差
C2 = diag([1/maxstd, 1/maxstd, 1]) #対角項
C2[0][2] = -m[0]/maxstd
C2[1][2] = -m[1]/maxstd
tp = dot(C2,tp) #内積

途中までは方法が同じですので説明を省き、最後の内積だけ開始点と同じように説明してみます。
dotは内積で下記が計算結果となります。
\begin{pmatrix}-0.6363961&0.77781746&0.77781746&-0.91923882\\-0.84852814&-1.13137085&1.13137085&0.84852814\\1&1&1&1\end{pmatrix}
f:id:Yaju3D:20140330214842j:plain
対応点の4点座標は左上から反時計周りで(20, 20)、(0,120)、(140,120)、(160, 0)となっています。
内積の計算結果で1行目がY座標、2行目がX座標とし対応点と位置を合わせてみます。
重心を原点0として対角線の長さで4点座標とすると下記のようになります。

(20, 20) (140,  0) → (-0.84852814,-0.6363961 ) (0.84852814,-0.91923882)
( 0,120) (160,120) → (-1.13137085, 0.77781746) (1.13137085, 0.77781746)
# 線形法のための行列を作る。対応ごとに2つの行になる。
nbr_correspondences = fp.shape[1]
A = zeros((2*nbr_correspondences,9))
for i in range(nbr_correspondences):
A[2*i] = [-fp[0][i],-fp[1][i],-1,0,0,0,
tp[0][i]*fp[0][i],tp[0][i]*fp[1][i],tp[0][i]]
A[2*i+1] = [0,0,0,-fp[0][i],-fp[1][i],-1,
tp[1][i]*fp[0][i],tp[1][i]*fp[1][i],tp[1][i]]

展開した行列

[[ 0.75        1.         -1.          0.          0.          0.   0.47729708  0.6363961  -0.6363961 ]
 [ 0.          0.          0.          0.75        1.         -1.   0.6363961   0.84852814 -0.84852814]
 [-0.75        1.         -1.          0.          0.          0.   0.58336309 -0.77781746  0.77781746]
 [ 0.          0.          0.         -0.75        1.         -1.  -0.84852814  1.13137085 -1.13137085]
 [-0.75       -1.         -1.          0.          0.          0.   0.58336309  0.77781746  0.77781746]
 [ 0.          0.          0.         -0.75       -1.         -1.   0.84852814  1.13137085  1.13137085]
 [ 0.75       -1.         -1.          0.          0.          0.   0.68942911 -0.91923882 -0.91923882]
 [ 0.          0.          0.          0.75       -1.         -1.  -0.6363961   0.84852814  0.84852814]]

特異値分解をします。

#特異値分解
U,S,V = linalg.svd(A)
H = V[8].reshape((3,3))

特異値分解したVの値が最小二乗法の解となります。

[[ 0.58620495 -0.04611386 -0.05921439]
 [-0.01117912  0.55336629 -0.05869036]
 [-0.10869141 -0.05928622  0.57062988]]

値を調整します。逆行列を使っているが、何をやっているか調査中

# 調整を元に戻す
H = dot(linalg.inv(C2),dot(H,C1))

調整した値です。

[[  4.29825103e-01  -8.89293317e-02   1.42286931e+01]
 [ -1.18572442e-01   4.29825103e-01   1.42286931e+01]
 [ -1.35864257e-03  -7.41077764e-04   7.11434654e-01]]

H[2,2]を1.0にするため、値を調整します。

# 正規化して返す
return H / H[2,2]

正規化した値です。

[[  6.04166667e-01  -1.25000000e-01   2.00000000e+01]
 [ -1.66666667e-01   6.04166667e-01   2.00000000e+01]
 [ -1.90972222e-03  -1.04166667e-03   1.00000000e+00]]