前記事で、「homography.pyでは、更に点群を平均値が0、標準偏差が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となります。
これは何をしているのか？
アフィン変換に拡大縮小(sx,sy)と平行移動(tx,ty)があります。
これに当てはめると「1/標準偏差の最大値」で拡大縮小させた後に「平均値/標準偏差の最大値」で平行移動させていると思われます。
この部分が本の記述にあった「平均値が0に、標準偏差が1になるように正規化」にあたるのではないでしょうか。
fp = dot(C1,fp) #内積
dotは内積で下記が計算結果となります。
※内積するのは開始点(fp)にアフィン変換行列(C1)を掛けた行列を求めるためです。
参考:文系のための「内積」（２）
この値が何を示しているのを下図で説明します。
この図は、縦幅(120)と横幅(160)から最大値である横幅(160)で正方形にして、赤線が交わったところが重心となります。ちなみに縦幅(120) / 横幅(160)=0.75となります。
開始点の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は内積で下記が計算結果となります。
対応点の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]]
