誤差を最小化する相似変換行列の求め方

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

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

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

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

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

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

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

No.4 回答者： stomachman 回答日時： 2014/09/11 14:00

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

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

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

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

No.3 回答者： stomachman 回答日時： 2014/09/10 10:01

ANo.2訂正。

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


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

に訂正です。

No.2 回答者： stomachman 回答日時： 2014/09/10 00:29

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

これはご質問にある通り、繰り返し数値計算をやって最適解を探索するしか手がありません。パラメータ（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回繰り返す」と決めとく、なんてのも、実用上はアリです。

この回答へのお礼

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

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

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

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

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

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

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

No.1 回答者： trytobe 回答日時： 2014/09/09 23:16

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

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

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

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

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

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

ポリゴン 法線ベクトル 外積 - 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 …

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

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

この回答へのお礼

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

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

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

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

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

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