三次元の変換行列について質問です。

空間上のある3点について、座標変換前と変換後の座標が与えられているとします。
その座標を元に、行った変換の表す行列を知りたいと考えております。
ただし変換は並進と回転のみとし、剪断変形等は行わないという条件です。
(行列は同次変換の形などで得られれば良いです)

ただし、変換前と変換後では三点の相対的な位置関係は完璧には一致しておりません。
そのため得られた変換を元の3点に行った際に、変換後の3点の座標(既知)との
誤差を最小化したいという問題です。

変換が相似変換でなければ疑似逆行列により計算が可能と思いますが、
並進・回転のみなのでどうしようか聞きたいと思い質問をしました。

パラメータは「x, y, z 軸方向の並進移動量 + x, y, z 軸周りの回転量」の6つなので、
プログラムを組んで6パラメータを少しずつ変化させ、
誤差の二乗和が最小となる箇所を探すという方法は思いついております。
ですが、その他に行列計算により解析的に答えを得る方法はないでしょうか。

また、変換前後の3点を元に誤差を最小化する変換行列を得ることが可能であるなら、
座標の情報が3点でなくそれ以上の場合でも計算は可能でしょうか。

よろしくお願いいたします。

このQ&Aに関連する最新のQ&A

A 回答 (4件)

ANo.2へのコメントについてです。



> 例えば評価を誤差の絶対値の和などにすれば解けるということはあるのでしょうか。

 三角関数が入っている以上、非線形性は取り除けません。それでも何か旨い手はないかと、かつていろいろ考えたこともありますけれども、結局は反復計算しかありませんでした。ならば、下手な工夫をやるよりも、品質の良い出発値を使って素直に反復する方が計算が速くできます。

 なおご質問は、元の点のそれぞれが、座標変換によってどの点に移ったのかが分かっている、という話ですから簡単な部類です。もし、移った先が一体どの点なのか、変換前後の点同士の対応は不明、ということだと、なかなか手のこんだことをやらねばなりません。
    • good
    • 0

ANo.2訂正。



> (たとえばy[1]-y[2])が変換で(z[1]-z[2]に)一致するように回転を決める。


(たとえばy[1]-y[2])が変換で(z[1]-z[2]に)平行になるように回転を決める。

に訂正です。
    • good
    • 0

 誤差の二乗和を評価に用いていますから、非線形最小二乗法の問題です。

これはご質問にある通り、繰り返し数値計算をやって最適解を探索するしか手がありません。パラメータ(6次元ベクトルxで表しましょう)で表される変換P(x)の推定値を、繰り返しのたびに少しずつ改良して行く訳です。(変換Pを指して推定「値」と呼ぶのは日本語としてちょっと変な感じがしますけれども、あらゆる変換がなす空間の中での1点がPだ、というキモチを表しています。)

 解析の対象になるのは、繰り返しごとに「変換の現在の推定値」P(x) から「変換の改良された推定値」P(x+Δx) をどうやって推定すると早く正しく収束するか、また簡単に計算できるか、ということだけです。
 どうやれば早く(少ない繰り返し回数で)正しく収束するか、ということについては、非線形最小二乗法の一般論として様々な工夫が開発されています。しかし、ご質問の状況では、問題の性質が素直であり、しかも計算量がごく小さいので、そんな工夫は必要ない。むしろ凝ったことをやって毎回の繰り返し計算が複雑になることの弊害の方が大きいでしょう。

 以下で提案する方法は「ガウス・ニュートン法」の一種で、概ね「繰り返しのたびにパラメータの有効数字の個数(つまり有効桁数)が一定量ずつ増える」という性質を持っています。実用上、せいぜい数回~十数回の繰り返しをすれば足りるでしょう。

[1]「変換の改良された推定値」P(x+Δx) の計算には、線形近似を使います。既に得られている「変換の現在の推定値」 P(x)に、さらに微小な変換ΔPを掛けることで「変換の改良された推定値」P(x+Δx)を得るものとします。
  P(x+Δx) = (ΔP) P(x)
そして、その微小な変換ΔPをパラメータの一次式で近似するのです。すなわちP(x+Δx)をΔxでテイラー展開したときの二次以上の項を全部無視することで、
  ΔP ≒ Σ((∂P/∂x[i])(x))Δx[i]  (Σはi=1~6の総和。xは「パラメータの現在の推定値」)
とやる。ご質問の場合ですと具体的には、「ひとつの軸について微小な角度Δθだけ回転する」という操作を
  sin(Δθ) ≒ Δθ
  cos(Δθ) ≒ 1
で近似してやれば良いだけです。3つの回転それぞれについてこの近似を行うと、微小な変換ΔPは回転角度(3つ)と並進(3つ)のパラメータの線形結合(一次式)で表されたことになります。
 この線形近似によって、ΔPの6つの(微小な)パラメータΔx[i](i=1~6)の最適値は線形最小二乗法を使ってイッパツで計算できる形になります。すなわち変換前の点をy[k]、変換後の点をz[k]、残差をε[k] (k = 1, 2,…, K)として
  minimize Σ(|ε[k]|^2) (Σは全ての点 k = 1, 2,…, Kについての総和)where,
  z[k] = (ΔP) (P(x) y[k]) + ε[k]
という線形最小二乗法の問題を解くだけです。こうしてΔxが得られるので、「パラメータの改良された推定値」(x+Δx)が決まり、これを使って「変換の改良された推定値」P(x+Δx)を算出します。

[2] 実用上の非常に大切なポイントは、パラメータと変換の最初の推定値(出発値)をどう決めるかです。正解に近いところから出発すると、少ない繰り返しで安定して答が出るからです。
 ご質問の状況ですと、たとえば「3点が張る平面の法線ベクトルが変換で一致し、かつ、あるひとつの点から別のひとつの点へのベクトル(たとえばy[1]-y[2])が変換で(z[1]-z[2]に)一致するように回転を決める。さらに、点群の重心が変換で一致するように並進を決める」という問題を解析的に解いて、答を出す。この計算によって、「ほとんど正しい」パラメータの推定値および変換の推定値が得られますから、これを出発値として使い、繰り返し計算で改良していく訳です。

[3] どこで繰り返しをやめるか、ということも考えておかなくちゃいけません。前の繰り返しに比べて残差二乗和の相対変化がうんと小さくなったら、とか、パラメータの変化がうんと小さくなったら、など、用途に応じて適切な打ち切り条件を決めておく必要があります。が、その判定のための計算があんまり複雑じゃ本末転倒。逆に、大雑把に「ナニハトモアレ10回繰り返す」と決めとく、なんてのも、実用上はアリです。
    • good
    • 0
この回答へのお礼

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

「誤差の二乗和を評価に用いていますから、非線形最小二乗法の問題です」
とのことですが、例えば評価を誤差の絶対値の和などにすれば解けるということはあるのでしょうか。

誤差の二乗和というのは「できるだけ両者を一致させる」という目的のための
指標の一つですので、問題が解けるのであれば他のものでも問題ありません。
もしご存知でしたら、お教えいただければと思います。

ただ、三次元の回転行列は例えば
・sin(Alpha)・cos(Beta)・cos(Gamma) + cos(Alpha)・sin(Gamma)
のような成分が入るので、どのみち無理なのかなと思いました。

数値計算に関してもお教えいただき、ありがとうございます。
6つもあるパラメータをどう変化させて探索するか、
全探索だと時間がかかるだろうかなどと考えていたので、
このあたりを教えていただき助かります。

ガウス・ニュートン法についても調べてみて、
まずはそれらしいプログラムを書いてみようと思います。

お礼日時:2014/09/10 21:54

(こちらを情報工学カテゴリーに書いたほうが実践的だったかもしれませんが、)



変換前の3点を△ABC とみなし、変換後の3点を△A'B'C' とみなして、

それぞれの三角形の法線ベクトルを求め、

2つの法線ベクトルの向きが平行になるように法線回りの回転(ひねり角度)と、法線を飛行機の機首のように持ち上げる回転とを行い、

その後は、法線ベクトルが同じ点に揃うように平行移動(並進)させる、

というだけで十分ではないでしょうか。

ポリゴン 法線ベクトル 外積 - Google 検索
http://www.google.co.jp/search?q=%E3%83%9D%E3%83 …

法線ベクトル 外積 座標変換 - Google 検索
http://www.google.co.jp/search?q=%E6%B3%95%E7%B7 …

並進分は単純な引き算ですので、法線ベクトル(外積)の向きを平行に合わせる(法線ベクトル2つの「内積」を0にする)ような「ひねり回転」と「仰角・俯角」の2パラメータだけが試行錯誤(最小二乗法)になるだけで済むのではないでしょうか。

そして、3角形→3角形に限らず、4角形以上のポリゴンでも、同じ長さ・角度を維持している3点を用いて、その三角形の並進・回転で、ポリゴン全体の並進・回転を成立させることも可能かと思います。
    • good
    • 0
この回答へのお礼

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

私も一番最初に、三角形の法線をそろえてから計算するというのを考えました。
ですが、回答中の「法線ベクトルが同じ点に揃うように平行移動(並進)させる」
というのがいまいち分かりませんでした。

三角形をABCとA'B'C'とした時、確かに三角形の法線を求め、
それを一致させる回転を求めることはできます。
しかし、法線ベクトルはあくまでABとACの外積として求めるものであり、
「ベクトルの始点をそろえる」ということに意味はないのではないでしょうか。

できるのは、例えば点Aと点A'を一致させるとか、
ABCとA'B'C'の重心をそろえることだけではないかと思います。

また2つの三角形の法線をそろえる(二つの三角を同一平面内におさめる)のが、
もっとも最適なそろえ方になるのか、という点に自信がありません。
もしかしたら、それが最適なのかもしれませんが...

お礼日時:2014/09/10 21:32

このQ&Aに関連する人気のQ&A

お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!

このQ&Aを見た人が検索しているワード

このQ&Aと関連する良く見られている質問

Q回転した座標系を基準とし、再回転したときの回転行列について

x軸、y軸、z軸が互いに直角に交わる座標系を考えます。(これを座標系Aとします)
座標系Aを、原点を中心とし、各軸ごとにθxa,θya,θzaだけ回転させた座標系を座標系Bとします。
さらに、座標系Bを基準とし、各軸ごとにθxb,θyb,θzbだけ回転させた座標系を座標系Cとします。

このとき、座標系Aから見た座標系Cの回転角は、どのように計算すればよろしいでしょうか?

座標系Aを基準とした回転角で座標系Bを計算し、さらに座標系Aを基準とした回転角で座標系Cを計算し……という問題であれば、単純に回転行列を掛けていけばいいと思うのですが、
「1つ前の座標を基準とした回転角を与えられたとき、全体でどれだけ回転したか?」
を表現する方法がわからなかったので、ご教示いただければ幸いです。

何卒よろしくお願いいたします。

Aベストアンサー

[1] 回転を組み合わせることについて

> 座標系Aを、原点を中心とし、各軸ごとにθxa,θya,θzaだけ回転させた座標系を座標系Bとします。

 ご質問では、どうも、これを一度の回転とお考えのように見受けられます。(違ったら失礼。)
 しかし、正しくはそうではない。x軸のまわり、y軸のまわり、z軸のまわりと3回の回転を組み合わせたんです。つまり、「さらに、座標系Bを基準とし、…」を持ち出すまでもなく、もうすでに、複数回の回転を組み合わせたものをお考えなのです。
 そのうえ、この文章だけではどんな回転をしたのかさっぱり分からない。というのは、第一に、いろんな回転を何度も繰り返す場合、(ご承知の通り)やる順番を変えれば結果も変わるからです。

(1) これら三回の回転はこの順番でやったのかどうか。

 ま、仮にこの順番でやったのだとしましょう。で、最初にやったx軸のまわりでの角度θxaの回転は良いとして、次にやった回転は、

(2)「座標系Aのy軸」のまわりでの回転なのか、それとも、「最初にやった(x軸のまわりでの角度θxaの)回転の結果得られた座標系のy軸」のまわりでの回転なのか。

 これがはっきりしません。三番目の回転についても同様です。

[2] 回転の表現について
 上記[1]の曖昧さについては補足を求めません。なぜなら、「座標系Aから座標系B、および座標系Bから座標系Cへの変換を(曖昧な文章ではなく)行列で表現したらどうなるか」について、質問者は先刻ご承知のようだからです。では、そのご承知の内容を確認しましょう。
 原点を通る直線を中心軸とする回転は、関係

 R' R = I  (「'」は転置)

を満たす3x3の行列(直交行列)で表現されることはご存知の通り。逆に、この関係を満たす行列Rは、「原点通る直線を中心軸とする回転を行う操作」か、あるいは、「鏡に映すように反転してから原点通る直線を中心軸とする回転を行う操作」かのどちらかを表しています。

 さて、このような回転を幾つ組み合わせようとも、

> 単純に回転行列を掛けていけばいい

と仰る通りで、「直交行列R1,R2,…,Rnで表される回転を、この順番で適用する操作」をRとすると、
 R = Rn … R2 R1
となる。確かに
 R' R = (Rn … R2 R1)' Rn … R2 R1 = R1' R2'…Rn' Rn … R2 R1 = I
を満たします。
 ところで、Rは、「原点通る直線を中心軸とする回転」ですから、その直線の方向を表す単位ベクトルxがある。つまり、Rは単位ベクトルxのまわりでの回転を表しているわけです。回転変換を表す行列Rを与えた時、このベクトルxは回転によって変化しないのだから、
 R x = x
を満たします。このxをRの固有ベクトルと言います。
 Rを与えた時にxを知るにはこの方程式を解けば良い。これで(座標系Aにおける)回転軸の向きが分かります。一方、角度θは何かというと、「xと直交する適当なベクトルvと、それが回転Rによって移ったものRvのなす角度」のことですから、両者の内積を取って
  v' R v = cosθ
から計算できます。
 逆に、xとθが与えられたときにRを構成するには、「直交行列Rであって、R x = x を満たし、かつ、x' v = 0 となるような単位ベクトルvについて v' R v = cosθ を満たすもの」を考えればよい。

 (うるさいことを言うと、回転の中心軸の方向を表すベクトルは当然2つある。つまり互いに逆向きの単位ベクトルです。一方、回転角についても、どっちまわりをプラスとみなすか、のやりかたが2通りある。ですが、ま、そういう細かいことは教科書に任せます。)

[3] ご質問に戻って

> 「1つ前の座標を基準とした回転角を与えられたとき、全体でどれだけ回転したか?」
> を表現する方法

を文字通り(とは言っても不足の部分は補って)解釈すれば、「あるベクトルxと、そのまわりで回転した角度θを与えた時、xとθは?」という問いに他なりませんから、答は初めからそこにある。これじゃ質問になってない訳です。

 一方、(おそらく)ご質問の意図は、Rを「各軸ごとにθxa,θya,θzaだけ回転させ」る、という形式で表現したいということなのでしょう。そういうことを考えるためには、まず[1]で申し上げた曖昧さをきちんと整理する必要がある。その上で、Rを三つの回転の積で表すことを考えれば良い。

 しかし、そんな面倒な表現を使わねばならない場合は滅多にない。単にR、もしくはxとθで表した方が単純明快だからです。

[1] 回転を組み合わせることについて

> 座標系Aを、原点を中心とし、各軸ごとにθxa,θya,θzaだけ回転させた座標系を座標系Bとします。

 ご質問では、どうも、これを一度の回転とお考えのように見受けられます。(違ったら失礼。)
 しかし、正しくはそうではない。x軸のまわり、y軸のまわり、z軸のまわりと3回の回転を組み合わせたんです。つまり、「さらに、座標系Bを基準とし、…」を持ち出すまでもなく、もうすでに、複数回の回転を組み合わせたものをお考えなのです。
 そのうえ、この文章だけでは...続きを読む

Q行列 同次変換 回転 並進

行列 同次変換 回転 並進

同次変換について質問させて下さい。
同次変換はよく、(X Y Z 1)のように
4列の行列で表されます。

4行4列の同次変換を表す行列で、
例えば、4行目が(1 2 3 1)とはXに1、
Yに2、Zに3の並進を表しています。


ここで質問なのですが、最後の4列目
はなぜ1と書かれるのでしょうか?

Aベストアンサー

1 でなくてもよいのですが、1 だと便利なので、そうしています。
4数対 (x y z w) によって、空間の点 (x/w, y/w, z/w) を表す
のが同次座標です。この表示法だと、回転、鏡影、拡大縮小の他に
並進まで行列の乗法で表すことができ、計算に便利なのですが、
(x y z w) を位置ベクトルに翻訳するのに、いちいち /w の除算を
しなければならないのは面倒です。w = 1 に固定してあれば、
そのような手間がなくて楽なのです。

Q「行列Aの固有値=1」⇒「Aは回転変換を示す行列」?

以前こちらでお世話になりましたaonegiです。
先日タイトルどおりの質問をさせていただき回答をいただいたのですが、よく吟味すると真に知りたいことが求められなかったのでまた新たに質問を出させていただきました。(前回の内容⇒http://oshiete1.goo.ne.jp/kotaeru.php3?q=1733069)

早速質問に入らせていただきます。
問題なのは行列Aの固有値をλとした時、|λE-A|=0を計算して固有多項式が求まるハズなんですが、そこで固有多項式の解がλ=1の場合、行列Aが回転を示す行列らしいのですが、何故そう言えるのか分かりません。
申し訳ありませんが再度御教授願います。

Aベストアンサー

はっきり言えば、間違っています。
回転を表す行列以外で固有値が1の行列はたくさんあります。
また、回転を表す行列の固有値の絶対値は1ですが、固有値そのものは1ではないことがほとんどです。
二次元では恒等変換以外の回転変換の固有値は1ではありません。
三次元では回転軸方向の固有ベクトルに対する固有値は1なので固有値1を確かに持っています。(1以外の固有値も持っているが)
一般的に対角線に1を並べてそれより下側はすべて0、上側には好きな数字を入れると固有値1の行列ができます。しかしながら、特殊な場合以外は回転行列になりません。
何かを勘違いされていると思います。

Q座標系 三角関数 回転行列

座標系 三角関数 回転行列

単位円において、θの正方向は反時計回りですがこれって右手系の座標系を採用しているからですか?

x,y,zの3軸を考えると以下になってしまって混乱しています。



1.(Z軸が奥向きの場合)
Y



 --→X
この座標系では、時計回りがθの正の方向ですよね?



2.(Z軸が手前向ききの場合)
   Y
   ↑
   |
   |
X←--
この座標系では、反時計回りがθの正の方向ですよね?

参考書などに書かれている単位円はθの正方向は反時計回りだと認識しているのですが、その場合、2の座標系なんでしょうか?

私は、1の座標系で反時計回りを正にしているように思いました・・・

回転行列を作ろうと思って考えた結果、自分でも混乱しています。。。

申し訳ないのですが、ご回答よろしくお願い致します。

Aベストアンサー

#2,#5,#6です。

A#6の補足の質問について回答

>つまり、右手系でX→Y→Z→Xと輪環の順に重ねる方向を正と定義すれば、XY平面、YZ平面、ZX平面で軸は全て原点手前方向になりますね。

その通り。

>こうすると、全て表面から見ていることになり反時計回りがθの正方向と一貫しています。
この認識で良いでしょうか?

その通り。その認識で結構です。

これで輪環の順に座標軸を重ねるように座標軸を回転させたとき、残りの座標軸が上(手前)の方を向いていれば反時計回りがθの正方向、輪環の順に座標軸を重ねるように座標軸を回転させたとき残りの座標軸の正方向が下(裏側)の方を向いていれば、θの正方向は時計回りになる(反対側の裏から見れば反時計回りになる)というわけです。

こうして一貫して覚えていれば、θの正方向を間違えることはなくなるでしょう!

Q座標上のある点が、ある3つの座標点で結んだ三角形の領域内にあるか調べる

座標上のある点が、ある3つの座標点で結んだ三角形の領域内にあるか調べる方法。

座標上に3つの点(x1,y1)(x2,y2)(x3,y3)で結ばれた三角形があります。
ある点(px,py)が、この三角形の内側の領域に存在するかどうかを知りたいのですが、
数学のなんという分野で、どういう求め方をするのかがわかりません。
どなたかお力添えいただければ幸いです。
関係ないかもしれませんが、左上を0,0とし、右下はn1,n2の、
Windowsペイントのようなマイナスを考慮しない座標になっています。
線上を内側とするか、外側とするかはどちらでもかまいません。
どなたかお詳しい方、お暇なときにでもご回答よろしくお願いします。

Aベストアンサー

点(x2,y2)を視点にして点(x3,y3)の方向を見たとき、点(xp,yp)が右側にあるか、左側にあるか
点(x3,y3)を視点にして点(x1,y1)の方向を見たとき、点(xp,yp)が右側にあるか、左側にあるか
点(x1,y1)を視点にして点(x2,y2)の方向を見たとき、点(xp,yp)が右側にあるか、左側にあるか
この3つの位置する方向がすべて同じなら、点(xp,yp)は3点で結ばれた三角形の内側にあります。

式で表すと、
z1=(x3-x2)(yp-y2)-(y3-y2)(xp-x2)
z2=(x1-x3)(yp-y3)-(y1-y3)(xp-x3)
z3=(x2-x1)(yp-y1)-(y2-y1)(xp-x1)
を計算して(外積の応用)、
z1>0,z2>0,z3>0 または z1<0,z2<0,z3<0 なら三角形の内側
z1=0,z2>0,z3>0 または z1=0,z2<0,z3<0 なら線分(x2,y2)-(x3,y3)上
z1>0,z2=0,z3>0 または z1<0,z2=0,z3<0 なら線分(x3,y3)-(x1,y1)上
z1>0,z2>0,z3=0 または z1<0,z2<0,z3=0 なら線分(x1,y1)-(x2,y2)上
z2=0,z3=0 なら点(x1,y1)上
z1=0,z3=0 なら点(x2,y2)上
z1=0,z2=0 なら点(x3,y3)上
上記以外は三角形の外側

点(x2,y2)を視点にして点(x3,y3)の方向を見たとき、点(xp,yp)が右側にあるか、左側にあるか
点(x3,y3)を視点にして点(x1,y1)の方向を見たとき、点(xp,yp)が右側にあるか、左側にあるか
点(x1,y1)を視点にして点(x2,y2)の方向を見たとき、点(xp,yp)が右側にあるか、左側にあるか
この3つの位置する方向がすべて同じなら、点(xp,yp)は3点で結ばれた三角形の内側にあります。

式で表すと、
z1=(x3-x2)(yp-y2)-(y3-y2)(xp-x2)
z2=(x1-x3)(yp-y3)-(y1-y3)(xp-x3)
z3=(x2-x1)(yp-y1)-(y2-y1)(xp-x1)
を計算して(外積...続きを読む


人気Q&Aランキング

おすすめ情報