円柱の抽出

3Dの点データから円柱をフィットして求めたい(最小二乗近似など)のですが、良い方法が見当たりません。
どなたか教えて下さい。
適当な参考書やURLをご存知ならば
合わせて教えて頂けるとありがたいです。

No.3 回答者： stomachman 回答日時： 2000/12/15 16:42

思いっきり非線形ですねえ。

３次元空間で
　　z軸を中心軸とする半径rの円筒を(自由度1)
　　平行移動し(自由度2)　--- z軸方向の移動は無意味なので1自由度減るから
　　回転し(自由度2)　　　--- 円筒の自転は無意味なので1自由度減るから
たものということですから、まっとうに行けば、5つの自由度(パラメータ)を決定する必要があります。r^2を線形最小二乗で決めるとしても、なお4次元空間の探索が必要です。それは最小二乗法の教科書を見ていただくことにして(「最小二乗法による実験データ解析」は名著です）、
正攻法の嫌いなstomachmanとしては、この問題を2次元空間の非線形探索に帰着します。

［１］仮に点列の回転(自由度２)ができ、この状態で点列
(Xn,Yn,Zn) (n=1,2,....,N) は
(Xn-Xc)^2 + (Yn-Yc)^2 = r^2 + εn ... (1)
という状態になったとする。（つまり、円筒の中心線がちょうどZ軸と平行になっている訳です。）r^2が円筒の半径(の二乗)の自由度、Xc,Ycが平行移動の自由度で、εnは誤差です。

［２］普通のやり方だとXc, Yc,r^2を決める手続きは非線形ですが、「誤差を小さくできるようなXc,Yc,r^2が存在する」と仮定すると、線形の問題に変換できる。すなわち誤差=0とおいた(1)式を整理すると
(2Xn)Xc+(2Yn)Yc +(r^2-Xc^2-Yc^2) = (Xn^2 +Yn^2)
ここで
Xn,YnおよびWn=(Xn^2 +Yn^2)は与えられたデータである。
Xc,Yc,およびK = (r^2-Xc^2-Yc^2)が未知数である。
これで線形最小二乗法の問題になりましたので、一発で解けます。すなわち、行列で書くと
A p = q
ここに既知のN行3列の行列AのA[n,m] (n行m列成分)はA[n,1] = 2Xn, A[n,2]=2Yn, A[n,3] = 1
未知の３次元縦ベクトルpのp[k](k行成分)は p[1]=Xc, p[2]=Yc, p[3]=K
また既知のn次元縦ベクトルqのq[n](n行成分)は q[n]=Xn^2 +Yn^2
計算速度や精度をあんまり気にしないでこれを解くと
p=(A' A)~ A' q
ここに A'はAの転置行列、(A' A)~は(A' A)の逆行列です。(この計算はExcelでも出来ますよ。関数minverse, mmult, transposeを使います。）
Xc,Yc,Kからr^2を出すのは簡単ですね。（知りたいのはrではなくr^2）

［３］さて、こうして決めたXc,Yc,r^2を使って、(1)式の誤差の二乗和を計算します。
S = |ε1|^2 + |ε2|^2 + ..... +|εn|^2
もし［１］で行った回転が旨くいっていれば、Sは非常に小さくなるし、さもなければSは巨大になる。
だから、Sが最小になるような回転を探せば良いのです。

［４］さて、回転です。
元の点列 (xn,yn,zn) (n=1,2,....,N)を3行N列の行列P
P[1,n]=xn, P[2,n]=yn, P[3,n]=zn
で表したとき、回転された点列は(Zは不要なので）2行N列の行列Q
Q[1,n] = Xn, Q[2,n] = Yn
で表されます。
x軸まわりの回転を行列 Rx=
１　０　０
０　Cx　-Sx
０　Sx　Cx
によって行い、引き続きy軸まわりの回転、およびZ成分の削除を行列 Ry=
Cy　０　-Sy
０　１　０
でやるとすると
Q = Ry Rx P
つまり(Ry Rx) =
Cy　-SxSy　-CxSy
0　　Cx　　　-Sx
ですね。ただし
Cx = cos(α), Sx= sin(α）
Cy = cos(β), Sy=sin(β)
ですから２個のパラメータα、βがある。この２つの数値をいろいろ変えて、Sを小さくするものを探索することになります。高速にやるため、あるいは安定化するためのいろんなテクニックが知られていますが、教科書を読んで貰うとして、ここでは手抜きしましょう。
[step1] (α,β)をとにかく決めて、Sを出してみる。これをS0とする。
[step2] Δα,Δβを適当な小さい値とする。
[step3] (α+Δα,β)でSを出してみる。
S0より大きくなったら、Δαの符号を逆にして、Sを出してみる。
それでダメならΔαを半分にしてやりなおす。
それでダメならΔαをさらに半分にして.....
とにかくS0より小さいSが出るまでやる。これをS1とする。この時のΔαを憶えておく。
今度は(α+Δα,β+Δβ)でSを出してみる。
S1より大きくなったら符号を変え、それでダメならΔβを半分にしてやりなおす。...
とにかくS1より小さいSが出るまでやる。この時のΔβを憶えておく。
これで、初めの(α,β)よりも良い(α+Δα, β+Δβ)が出た。
[step4] step3で最終的に得たΔα, Δβの3倍をΔα, Δβの値とする。そしてstep3へ。
（半分とか3倍とか、はまあ、いい加減な値で良いんです。）

［５］出発値の選び方
(α,β)=(0,0)とか90度とか、あんまりチョッキリの値はやめましょう。特異姿勢というのに嵌ってしまうおそれがある。
CG用レンダリングソフトでも使って、だいたいの値が分かる、というのでもよい。
また、たとえば点が円筒上に概ね一様に散らばっているというのなら、だいたいの向きを主因子法で決められます。(P P')をスペクトル分解し、固有値と固有ベクトルを求めます。円筒の半径に比べて（点の存在する範囲の）長さが長い、ということが分かっているなら、最大の固有値に対応する固有ベクトルの向きが円筒の向きの近似値です。逆に円筒の長さは明らかに短いというのなら、一番小さい固有値の固有ベクトルを取ればよい。

●円筒のフィッティングを１回やるだけなのか、大量のPが与えられて、次々自動的に円筒を決めていく必要があるのか、それによっても、計算法の根性の入れどころが変わってきます。がんばって-

この回答へのお礼

非常に丁寧な解説ありがとうございます。
たすかります。
ゆっくり考えてみます。

お礼日時：2000/12/15 18:49

No.2 回答者： stomachman 回答日時： 2000/12/14 18:40

(1)どんな円柱が来るのか分からないのか

(2)円柱の軸の向き（３次元的方向）があらかじめ分かっているか
(3) (2)に加えて円柱の軸の位置がわかっているか（つまり円柱の半径だけ分かれば良いのか）

によってしんどさが随分違います。この区分を補足していただけませんか？

この回答への補足

点群のデータのみです。
つまり点群データに直に円柱をフィットさせたいのです。

補足日時：2000/12/15 10:41

No.1 回答者： tullio 回答日時： 2000/12/14 18:38

一般化ハフ変換か，超2次曲面の当てはめをしてはいかがでしょうか．どのような最適化手法を使うかが工夫のしどころになります．