３次元空間内での線分の交差判定について

質問

はじめまして。
３Ｄ関係のプログラムを組む上で、線分同士の判定を行う必要があるのですが
数学の知識が乏しく困っています。

３次元空間内の線分ABとCDが交差しているか判定し、
交差していればその交点を求めたいのです。
２次元の場合はできたのですが、３次元になるとどうやって計算すればよいのか
わかりません。（交差以外に、ねじれの位置関係があるんですよね？）
どなたか教えていただけると助かります。

回答No.2

ametsuchi · ベストアンサー

１）一般に「捩れの位置」になりますから、互いに最短の位置を求める問題に帰着します。a-kumaさんの言われるような連立一次方程式では未知数が２つ、式がx,y,z３つとなるので解けません。2次元ならyosizoさんの言われるように未知数と式の数が２つで簡単に解けます。

２）線分AB、CDの何れかの長さが<ε以下、または、線分AB、CDはほぼ平行ならゼロ割を起こすか答えが求まっても答えの精度が著しく低下します。このチェックも実際のプログラムでは絶対に必要です。

３）で、以下、２）のチェック済みであると仮定して.....。

４）たまたま最近３D－CAD用に作ったもので数式の求めた方を忘れてしまったが、

// t0:AB方向単位ベクトル、t1:CD方向単位ベクトル、として、
// float spd = <t0・t1>即ち、ベクトルt0とt1の内積

float det = spd*spd - 1.0f; // spd:AB方向単位ベクトル
v01 = C - A; //３Dベクトル (C-A)
u0 = (spd*<v01・t1> - <v01・t0>) / det;
u1 = (<v01・t0> - spd*<v01・t1>) / det;

として、パラメータ、u0,u1が求まります。ここに、
u0:AB方向パラメータ、u0=0の時、q0(u0)=Aで、u0はAからの距離を表わす。
■q0(u0)=A+u0*t0...............Equ.1)
u1:CD方向パラメータ、u1=0の時、q1(u1)=Cで、u1はCからの距離を表わす。
■q1(u1)=C+u1*t1...............Equ.2)

これで、捩れの位置に２点求まるのですが、

・２点が距離の許容誤差よりも離れていたら、エラーとする。
・２点が距離の許容誤差以内なら、２点の中点を取る。
・中点がいやなら、重みを掛ける（場合もある）。

この説明で不十分ならあとで補足します。

回答No.7

ametsuchi

motsuanさん、アドバイスありがとうございます。

>外積を直接とっちえばいいものを ..中略...のような周りくどいことをやらないと 答えがもとめられないのかちょっと不思議に思います

一般に外積計算は内積計算に比べて「コストのかかる」計算であり、なるべく使わないようにしているだけです。motsuanさんのやり方の分かり易さや、「外積計算」の有用性については百も承知しているつもりです。

もし、私のやり方の方が「コストのかかる」計算であったなら、折角今回初めて「専門家」を自称したのにおハズカシい限りで、訂正させてもらいます。

motsuanさんのやり方も私がいう「一般性のある」やり方で、敬意を表します。３D幾何計算ライブラリなどでは直線（線分）同士は、「交点」ではなく「最近点」として求めるのが普通です。

お礼

返事が遅くなり、失礼しました。
詳しい解説、ありがとうございます。
＞３D幾何計算ライブラリなどでは直線（線分）同士は、「交点」ではなく「最近点」として求めるのが普通です。
そうなんですか、しりませんでした。勉強になります。

今回の皆さんの回答を参考にしてプログラムに落とし込んでみます。

回答No.6

motsuan

もうすでに完全に説明されているようですが、
なんとなく私も前に計算してなるほどと
なんか関心してしまったことがあるので
（他の方からみると見苦しいとは思いますが）
書きこませてください。

ametsuchiさんのアルゴリズムは
多分以下のように求めれば
すべてベクトル演算で求められると思います。
（ねじれの位置の場合も線分の中に対応する点が
　あれば一番近づいている点をあたえると思います。）

ベクトルを[v]と表します。
２つ直線[ｐ]（u）、[ｑ]（v）をパラメータ表示で表します。
それぞれ定点を[p0]、[q0]、大きさが任意の方向ベクトルを[s]、[t]とすると
[ｐ](u)＝[p0]＋[s]u
[ｑ](v)＝[q0]+[t]v
と表せます。このとき
[ｐ](u)－[ｑ](v)＝０
となる　u、v　があればそれが交点を与えます。
３次元空間で交わる場合には、
当然どの平面に射影したとしても(どっちの方角から見ても）
交点を結んでいなくてはなりません。
したがって適当な平面へ射影した場合の交点を求めれば良いわけです。
これがkent-mildsのおっしゃっていることで
やっぱり一番現実的なのでしょうね。

でも、なんとなく、
一番近づく点の特殊な場合が交点であって欲しい
という気持ちを収めるために以下のように基準面を選びます。
すなわち、一番近づく点が射影したときに交点となる平面は
[s]、[t]のベクトルで張る平面ですから、この平面を基準に考えます。
（以下、内積を・で、外積を×で表します。）

[ｐ](u)と[ｑ](v)を[s]、[t]のベクトルで張る平面に射影すると
１点に重なるためには[ｐ](u)－[ｑ](v)が面の法線と平行でなくてはなりません。
したがって、([ｐ](u)－[ｑ](v))×[s]と([ｐ](u)－[ｑ](v))×[t]は
面に平行なベクトルになりますから、これが[s]×[t]と直行する必要があります。
すなわち、
([s]×[t])・(([ｐ](u)－[ｑ](v))×[s])＝０・・・（＊）
([s]×[t])・(([ｐ](u)－[ｑ](v))×[t])＝０
が求める点の条件で、これから u, v が決まります。たとえば u を求めると、
([ｐ](u)－[ｑ](v))×[t]
　＝([p0]－[q0])×[t] + [s]×[t]u ＝０
両辺に[s]×[t]を掛けて
([s]×[t])・(([p0]－[q0])×[t] )＋([s]×[t])・([s]×[t])u ＝０。
これを解いて
u ＝－([s]×[t])・(([p0]－[q0])×[t] )/([s]×[t])・([s]×[t])。
同様に（tとsを入れ替えて）、
v ＝([s]×[t])・(([p0]－[q0])×[s] )/([s]×[t])・([s]×[t])
となります。

わたしはなんで
([ｐ](u)－[ｑ](v))と[s]×[t]の外積を直接とっちえばいいものを
（＊）のような周りくどいことをやらないと
答えがもとめられないのかちょっと不思議に思います。
以上です。お邪魔しました。

お礼

返事が遅くなり、失礼しました。
詳しい解説、ありがとうございます。
今回の皆さんの回答を参考にしてプログラムに落とし込んでみます。

回答No.5

ametsuchi

すいません。もう２つ補足。

１）a-kumaさんらの方法でもできないことはないのですが、一般性がなく、式３つのうちどの２つを選択するか決めなくてはいけません。私の方法なら、場合分けも不要で答えは一意的に求まります。

２）線分上に求まるか否かのチェックはa-kumaさんの言われるとおりです。パラメータの取り方が違っていますが、本質的な差異ではありません。私が示したパラメータだと、

■0<= u0,u1 <=線分長

だとOKです。

・しかし実際のプログラムでは、

■-tol<= u0,u1 <=線分長+tol

のようにすることが多いです。ここに,tolとは距離の許容誤差です。

※何か不明な点があったら聞いて下さい。

回答No.4

ametsuchi

補足します。

２線分上の点はパラメータで表現すると、

■q0(u0)=A+u0*t0...............Equ.1)
■q1(u1)=C+u1*t1...............Equ.2)

ですが、

・ベクトルq1(u0)-q0(u1)と、ベクトルt0、が直交
・ベクトルq1(u0)-q0(u1)と、ベクトルt1、が直交

という条件を課すと、２つのパラメータ、u0,u1に関する２元一次連立方程式になります。これを解いたのが、先程の、

float det = spd*spd - 1.0f; // spd:AB方向単位ベクトル
v01 = C - A; //３Dベクトル (C-A)
u0 = (spd*<v01・t1> - <v01・t0>) / det;
u1 = (<v01・t0> - spd*<v01・t1>) / det;

です。

お礼

すみません、これは↓のNo3　Kent-mild　さんへのお礼です。
＃ここの使い方に慣れてなくて、意味のない書き込みをしてしまいました・・・

返事が遅れて大変失礼しました。
みなさんからの回答を参考にしてプログラムを組んでみようと思います。
ありがとうございました。

回答No.3

kent-milds

プログラムで簡単に組める方法で説明します。
２次元ではできてらっしゃるようなので、２次元の説明は省略します。
３次元空間が（x,y,z）座標で定義されている場合、
xy平面上の交点のz座標が同一なら交差、一致しなければねじれです。
交点の座標はすぐでるのはわかりますよね？

お礼

３D幾何計算ライブラリなどでは直線（線分）同士は、「交点」ではなく「最近点」として求めるのが普通です。

回答No.1

a-kuma

線分AB上の任意の点Pはベクトル表現で

p = t×a ＋ (1 - t)×b
　　但し 0≦t≦1

です(本当は、記号の上に矢印を書きたいところなんですが)。
同様に、線分CD上の任意の点Qは

q = s×c ＋ (1 - s)×d
　　但し 0≦s≦1

です。線分が交差しているということは、PとQが同一になるという
ことですから、

p = q

を満たす t と s が有るかどうか、ということですね。
3次元であれば、各要素毎に3本の式ができ、変数が二つですから
簡単ですね。

例えば、x と y について、t と s についての連立方程式を解いて、
それぞれ０～１に入っていなければ、交差してない。
もし、その範囲に入っていたとすれば、z の式に代入してみて
等式を満たすようであれば、交差する、という判定でしょうか。

お礼

返事が遅くなり、失礼しました。
詳しい解説、ありがとうございます。
今回の皆さんの回答を参考にしてプログラムに落とし込んでみます。