この連立方程式の解き方を教えて下さい

質問

いつもお世話になっております．

下図の式を解きたいのですが解き方が分からなくて困っております．

これは，下図に記載されている行列Aを求める問題です．
（式1）を変形すると（式2）と（式3）が得られ，さらに変形すると（式4）と（式5）が得られます．
既知の条件として，（xs, ys）と(Xw, Yw, Zw)の組が6組以上分かっているとします．
このときに，行列Aを求めたいのですが，どうすれば解けるのでしょうか？

（式4）と（式5）に既知の条件を与えて，Ax=bのようにして逆行列を使って解きたかったのですが，
今回の場合，b=0となってしまい，この方法では解くことができません．（ここでの，A, b, xは下図の記号とは関係ありません．）

（式4）と（式5）において，a～lを変数として，最小二乗法を使ったりするのでしょうか？

数学にお詳しい方，解き方を教えて下さい．
ヒントや参考サイトでも構いませんので，宜しくお願い致します．

最終的には，この式を解くためのプログラムを組むことを考えています．
それでは，どうぞ宜しくお願い致します．

回答No.10

178-tall

>　Zc の推定は不可能なので、不可解 (解けない) なのかも…。

ANo.7 さんのご指摘に従って、C の代わりに
　3*1 行列 [xs；ys；1 ] を 4 列並べた 3*4 行列 Co
を使って
　CoW^(-1) = Bo
を求めたあと、たとえば各要素を l (小エル) で normalize し、
　3*4 行列 Bo/l =
　　　　　[ ao　bo　co　do ；
　　　　　　eo　fo　go　fo ；
　　　　　　io　jo　ko　l ]
　　　　　　　　　　　　↑
　　　　　これは「小エル」に非ず。数字の「イチ」。
でも勘定しますか…。

回答No.9

178-tall

>与えられる値は， （Xw，Yw，Zw） と （Xc，Yc，Zc ） の組ではなく （Xw，Yw，Zw） と （xs，ys） の組です

だとすれば、Zc の推定は不可能なので、不可解 (解けない) なのかも…。

回答No.8

hrsmmhr

問題にあっているのかどうか分かりませんが

計算出求まるセンサー値のデータ（回転・平行）はなんらかの座標系に従った値になるはずなのですが
この問題では入力からは、どの原点とどういう計測系で測った値として求めるのかが分からないので
比しか出てこないのではと思います。

少なくとも２点位での計測系での値（カメラで写されている物の実測とかを換算した結果）を定数項として用いれば定数項が出るはずです。どういう実測値を用いるかは、この計算を何に用いるかということによると思います

これが、何かの課題で、実測もできないのなら、適当な値を使えば値は出ますが、それがカメラのときに
何の意味があるのかが分からなくなります。課題の内容から推測できれば一番いいと思いますが。

回答No.7

nag0720

＃２です。

＃３さん、＃６さんの方法が有効でしょう。

Zcが分からないと・・・とありますが、
＃２の始めに書いたように、a～lは一意に決まるのではなく、比が決まるだけです。
なので、Xc=xs、Yc=ys、Zc=1としても問題ありません。

B=CW^(-1)
からa～lが求まりますが、その定数倍もまた解になります。

回答No.6

178-tall

ANo.3 さんのアプローチが有望らしい。少しばかりパラフレーズを…。

　3*4 行列 B =
　　　　　[ a　b　c　d ；
　　　　　　e　f　g　f ；
　　　　　　i　j　k　l ]

　4*1 行列 [Xw；Yw；Zw；1 ] の一次独立な 4 組 を 4 列並べた 4*4 行列 W
　対応する 3*1 行列 [Xc；Yc；Zc ] を 4 列並べた 3*4 行列 C

などの形式にして (式 1) を拡げてみると、
　C = BW　　　…(6)
ですね。

(6) の右から W の逆行列 W^(-1) を掛けたら、B が求まりませんか？
　CW^(-1) = B

お礼

178-tall様

解説ありがとうございます．
おかげ，No.3 さんの言ってることは
（条件）d=h=l=0以外はなんとなく分かりました．

ただ，与えられる値は，
（Xw，Yw，Zw） と （Xc，Yc，Zc ） の組ではなく
（Xw，Yw，Zw） と （xs，ys） の組ですので，
（xs，ys）から（Xc，Yc，Zc ）に変換する必要があるかと思います．

(式2)よりXc=xs*Zc
(式3)よりYc=ys*Zc
であり，教えて頂いた解法の前にZcを求める必要があります．

Zcは(式1)から分かる
Zc=iXw+jYw+kZw+lで求まりますが，この時点でi,j,k,lは未知です．

何か抜けていることがあるのでしょうか？

回答よろしくお願い致します．

回答No.5

hrsmmhr

最後の条件式は思い違いです
解ける気がしたんですが…
No.2さんのやり方でいいと思います

お礼

hrsmmhr様

回答ありがとうございます．

そうでしたか．
また何か分かることがありましたら宜しくお願い致します．

回答No.4

nag0720

＞i,l,j,kが未知数なので，i,l,j,kでわるような操作はできないと思うのですが，

実際に割る計算をするのではなく、
Xw(a/l)+Yw(b/l)+Zw(c/l)+(d/l)-xsXw(i/l)-xsYw(j/l)-xsZw(k/l)-xs=0
とし、
(a/l)、(b/l)などを未知数とみなして解くということです。

お礼

nag0720様

回答ありがとうございます．

そういうことだったんですね．
nag0720様の方法で解ける気がしてきました

先の補足からl=tzといことになります．
なのでtz=0が成立してしまうと，カメラと撮影対象が同一平面上に存在することになる気がするので・・・

先の補足から何か分かることがあれば，追加で教えて頂けると幸いです．
宜しくお願い致します．

補足

nag0720様

教えて頂いた方法ですと，lはどうやって計算するのでしょうか？
（xs, ys）と(Xw, Yw, Zw)の組が6組あれば，式を12本立てることができ，
このうち11本使い，教えて頂いた方法を使えば，AX=bの形で，
l≠0と仮定した場合に
a/l,b/l,・・・・,k/lは求めることができますが，lを計算することができません．
また，余った式からでは，上手く，変形することができません．

回答よろしくお願い致します．

回答No.3

muturajcp

(Xw,Yw,Zw)の組の中の一次独立な3組を
(xw1,yw1,zw1)
(xw2,yw2,zw2)
(xw3,yw3,zw3)
としてそれぞれに対応する(xs,ys)の組を
(xs1,ys1)
(xs2,ys2)
(xs3,ys3)

W=
(xw1,xw2,xw3)
(yw1,yw2,yw3)
(zw1,zw2,zw3)
C=
(xs1,xs2,xs3)
(ys1,ys2,ys3)
(1..,1..,1..)
とするとWの列ベクトルは一次独立だから
行列式|W|≠0だから
Wの逆行列W^{-1}が存在するからそれを求めて
B=CW^{-1}
を求めて、
B=
(a,b,c)
(e,f,g)
(i,j,k)

d=h=l=0
とすると

(xs1,xs2,xs3)=(a,b,c,d)(xw1,xw2,xw3)
(ys1,ys2,ys3).(e,f,g,h)(yw1,yw2,yw3)
(1..,1..,1..).(i,j,k,l)(zw1,zw2,zw3)
(1..,1..,1..).(0,0,0,1)(1..,1..,1..)

お礼

muturajcp様

回答ありがとうございます．

muturajcp様の説明のうち，
「Wの逆行列W^{-1}が存在するからそれを求めて」というところまでは分かりました．
それ以降が，私では理解することができません．
宜しければ，ご教授ください．

（１つ目）
B=CW^{-1}
B=
(a,b,c)
(e,f,g)
(i,j,k)
の２つの式はどういった条件で出したのでしょうか？

（2つ目）
d=h=l=0とすると
どういった条件で出したのでしょうか？

（3つ目）
(xs1,xs2,xs3)=(a,b,c,d)(xw1,xw2,xw3)
(ys1,ys2,ys3).(e,f,g,h)(yw1,yw2,yw3)
(1..,1..,1..).(i,j,k,l)(zw1,zw2,zw3)
(1..,1..,1..).(0,0,0,1)(1..,1..,1..)
の表記が何を示しているのかがわかりません．

回答よろしくお願いします．

回答No.2

nag0720

もしa～lが定数倍になったとしたら、Xc,Yc,Zcも定数倍になるだけだから、xs,ysの値は変化しません。
つまり、(xs, ys)と(Xw, Yw, Zw)だけではa～lは一意には決まりません。
決まるのは比です。

また、xs,ysが存在するということは、Zc≠0、よってi,j,k,lのうちどれかが0ではありません。
例えば、l≠0とした場合、
(式4)(式5)をlで割れば、Ax=bの形式にできるので解くことができます。

ただし、i,j,k,lのうちどれが0でないかは、与えられる条件によって変わってきます。
プログラムを組むとしたら、l≠0でうまくいかなかったら、k≠0の場合、j≠0の場合、i≠0の場合、とすべての場合を考慮しなければならないでしょう。

お礼

nag0720様

回答ありがとうございます．

xs,ysが存在するので，Zc≠0が成り立つということには気が付きませんでした．
たしかに，そうです．

ただ，今のところ，i,l,j,kが未知数なので，i,l,j,kでわるような操作はできないと思うのですが，
できるのでしょうか？

補足

muturajcp様，nag0720様，hrsmmhr様，また，その他の閲覧者の皆様

色々と考えて頂きありがとうございます．

問題を補足します．

この問題は，
世界座標系中の点(Xw,Yw,Zw)をカメラで撮影し，それの画像平面上での点を(xs,ys)としたときに，
その間に成り立つ変換式を求める問題です．

下図のAについて，もう少し詳しく説明すると
A=C[R|T]と表現できます．
Cはカメラの内部パラメータと呼ばれるもので
私が調べた限りでは，
C=
(fmx, s , mxPx, 0)
(0 , fmy, myPy, 0)
(0 , 0 , 1, 0)
(0 , 0 , 1, 1)
となり，この行列の要素は全て既知です．
（mx,my,Px,Pyについているx,yは添え字です．）

また，[R|T]はカメラの外部パラメータと呼ばれるものであり，
Rは回転行列，Tは並行行列であり，世界座標系からカメラ座標系に変換するための行列です
（r1,r2,r3,tx）
（r4,r5,r6,ty）
（r7,r8,r9,tz）
（ 0, 0, 0, 1）
と表記することができ，この行列の要素は全て未知です．
（上記同様にrについた数字とｔについたx,yは添え字です．）

結論をいうと，
この問題は，(Xw,Yw,Zw)と(xs,ys)の組が複数既知である場合に[R|T]を求める問題となります．

補足は以上です．
引き続きご教授して下さいますよう宜しくお願い致します．

回答No.1

hrsmmhr

例えば全部の条件式をlで割って、11次のAx=ｂの逆行列を計算して、
最後の条件式でlを決めるように計算できないでしょうか？

お礼

hrsmmhr様

回答ありがとうございます．
>最後の条件式
というのは，lでとりあえず割ってみて出てきた式ということでしょうか？
たしかに，lで割ることができれば，逆行列の計算で解くことができます．
ただ・・・，いまのところlが未知数です．

他にも条件が無いかを探してみます．
ありがとうございました．