「新・写真から立体を再現」では、エピポーラ方程式というものを使います（これについては別の回に説明します）。2枚の写真上の対応する点の組を方程式に代入してエピポーラ方程式に含まれる基礎行列Ｆを決定し、カメラ行列を求めます。このときＡｘ＝0の形の方程式の非自明解を求める必要があります。今回はこの非自明解を求める練習をしてみます。前に少しお見せしましたが、もう少し丁寧に。
基礎行列はランクが2と分かっています。そこでランク2の行列Ｄを用意し、Ｄｘ＝0の非自明解を求めます。Ｄｘ＝0 ⇔ Ｄ'Dｘ＝0（Ｄ'はＤの転置行列）なので、ここではＤ'Ｄｘ＝0の非自明解を求めることにします。まずＡ＝D'Dを計算し、Ａの固有値と固有ベクトルを求めます。固有値0に対応する固有ベクトルｖはＡｖ＝0ｖ＝0を満たすはずですから、固有ベクトルｖはＡｘ＝0の非自明解です。念のため、固有値に対応する固有ベクトルを求めた後Ａにかけて確かに0になることを確認します。さらに念のため、その固有ベクトルはＤｘ＝0を満たすことも確認します。数値計算なので、固有値や固有ベクトルの成分は十分小さければそれを0とみなします。
#========================================================
#固有値、固有ベクトル、非自明解
#
import sys
import numpy as np
D = np.array([[3, 1, 2], [1, 1, 1], [1, 1, 1]])
print('rank = ', np.linalg.matrix_rank(D))
print('')
A = np.dot(D.T, D) #A = DT*D
print('DT*D =\n', A)
print('')
B, C = np.linalg.eig(A) #固有値、固有ベクトルを求める
print('固有値 = \n', B)
print('')
#Cの列ベクトルが固有ベクトルなので、転置して取り出す
print('固有ベクトル = \n', C.T)
print('')
print('固有ベクトル = \n', C)
print('')
print('A * 固有ベクトル = ')
for i in range(0, 3):
print(np.dot(A, C.T[i])) #A*固有ベクトルを計算
print('')
print('D * 固有ベクトル = ')
for i in range(0, 3):
print(np.dot(D, C.T[i])) #D*固有ベクトルを計算
sys.exit() 実行結果は以下の通りです。
ーーーーーーーーーーーーーーーーー
rank = 2
DT*D =
[[11 5 8]
[ 5 3 4]
[ 8 4 6]]
固有値 =
[ 1.93808315e+01 6.19168480e-01 -1.97440383e-16]
固有ベクトル =
[[-0.74752958 -0.3638668 -0.55569819]
[-0.52395883 0.83721818 0.15662968]
[-0.40824829 -0.40824829 0.81649658]]
固有ベクトル =
[[-0.74752958 -0.52395883 -0.40824829]
[-0.3638668 0.83721818 -0.40824829]
[-0.55569819 0.15662968 0.81649658]]
A * 固有ベクトル =
[-14.48774491 -7.05204106 -10.76989298]
[-0.32441879 0.51837911 0.09698016]
[-8.8817842e-16 -8.8817842e-16 -8.8817842e-16]
D * 固有ベクトル =
[-3.71785192 -1.66709457 -1.66709457]
[-0.42139895 0.46988903 0.46988903]
[ 2.22044605e-16 -6.66133815e-16 -6.66133815e-16]
ーーーーーーーーーーーーーーーーー
うまく非自明解が求まっていることが分かります。
ここまでである程度Pythonで書くのにも慣れたし、以降は速いかも知れません。基準の立方体なしに立体の再現ができればこれで自分の中では一応ひと区切り、落ち着いてまた先へ進めます。繰り返しになりますが、こうして自分でコードを書くのとライブラリなどを使うのとでは達成感もまるで違うでしょう。もちろん目的によってはライブラリを利用した方がいいケースだってあると思いますが、理屈を勉強することが大切であることは変わらないと思います。
