ある形式で与えられる空間中の２直線間の距離・中点を求めたいのですが

質問

「空間中の２直線間に、共通垂線を一本引いて、垂線が各直線に対して持つ二つの足の間の距離(---1)と、二つの足の間の中点の位置(---2)」を、求めたいのですが、どのようにするのでしょうか.　

直線は、"平面xy座標上の一点aを通り、かつ、方位角bと仰角cとで定まる方向"という形式で与えられるのです.

最終目標としては、空間中にいろいろと与えられる二直線について、冒頭の「」を求めるプログラムを書きたいのです.

また、冒頭の「」内のことは、捩れの位置で始めて成立するわけですが、交わるのか捩れの位置なのかを、プログラムでもって、簡単に識別する方法ってあるのでしょうか?

このあたりの領域に詳しい方、どうぞ教えてください.

回答No.1

hogehogeninja · ベストアンサー

ベクトルを考えると考えやすいです。

まず、２直線をA,Bとし、それぞれの単位方向ベクトルをa^b^します。
（方位角と仰角からすぐにだせます）

また、xy座標平面上で、直線Aの交点 Aa から直線Bの交点 Ba に引いたベクトルをc_ 置きます。（ちょっと分かりにくいですが、x^ は単位ベクトル、x_ はベクトルを表すとします）
また垂線に平行な単位ベクトルをn^ と置き、２直線の距離（足の間の距離）をd と置きます。

すると
　n^⊥x^, n^⊥y^　　（１）
n^ = x^ × y^ （外積）　（２）
としよてく、また、α、βを直線A,B上で足の位置を表すパラメータとして
α x^ + d n^ = c_ + β y^　（３）
が成立します。この式の両辺似たいしてn^ との内積をとり、（１）とn^ の大きさは１であることに注意すると求める距離d は
　d = c_・n~　
と簡単に求まります。

また、（３）式からx^ と内積をとった式、y^ と内積をとった式、を作るなどしてαを求めます。

直線A のxy平面上の点Aa に引いた位置ベクトルをe_ で表せば、中点の位置は
　e_ + α x^ + (d / 2) n^
で表されます。

このやり方で求めたとき、交わるときには
　d = 0
を与え、特別にチェックする必要は無いです。
むしろ特別な配慮が必要なのは２直線の平行度で、平行な時には上の量は本当に数学的にも定義できなくなります。
平行に近いとき（厳密に平行でなくても）、プログラム上では非常に誤差が大きくなりますので、適当に誤差を見積もり、評価してやる必要があります。
（x^ と y^ の内積をとって平行度を調べるなどして）

お礼

ありがとうございます.最初はわからなかったのですが、解ったらすごく簡単に書けることに驚きました.以下のように、自分なりに、まとめてみたのですが、あっているでしょうか.
-------------------------------------------------------------------------------------
直線Aとxy平面との交点 a_=(ax, ay, az)
直線Aの方位角 θa
直線Aの仰角 φa

直線Bとxy平面との交点 b_=(bx, by, bz)
直線Bの方位角 θb
直線Bの仰角 φb

とすると、上記の値だけを使って、以下の左辺の値は、右辺のように書くことができる.

直線Aの単位ベクトル a^=(cos(θa)cos(φa), sin(θa)cos(φa), sin(φa))
直線Bの単位ベクトル b^=(cos(θb)cos(φb), sin(θb)cos(φb), sin(φb))

共通垂線の単位ベクトル n^=a^×b^

共通垂線の足間距離　d= (b_-a_)・(a^×b^)
[=即ち直線間距離]　　= (b_-a_)・n^

共通垂線足間中点位置　m_=a_ + α a^ + (d/2) n^
[即ち直線間中点]　　　　=a_ + α a^ + (1/2){(b_-a_)・(a^xb^)}・(a^×b^)

但し α=[{a^・(b_-a_)}-{b^・(b_-a_)*(a^・b^)}]/{1-(a^・b^)(a^・b^)} (但しa^・b^=1でない時)
-------------------------------------------------------------------------------------

回答No.5

ametsuchi

くどいですが、もう一つ補足。

２直線の外積を求めるやり方でも、私流のやり方でもコストは変わらんと言いましたが、実際に計算してみると、私流のやり方の方がコストが安いと分かりました。

直線上のパラメータt0,t1、v0,v1：直線の方向（単位ベクトル）、p0,p1：通過点とすると、直線上の点は

P0(t)=p0+t0*v0
P1(t)=p1+t1*v1

1) ２直線を通過する２点の差分ベクトル：p01=p0-p1
2) 内積 s01=v0・v1
3) 内積 s0 =v0・p01
4) 内積 s1 =v1・p01

として、
5) det＝s01^2 – 1.0
6) t0 = (-s0 + s01*s1) / det
7) t1 = ( s1 - s01*s0) / det

直線上のパラメータt0,t1が求まります。ここまでで、乗除算：11回、加減算：10回。

一方、外積を取るやり方だと、外積計算だけで、乗除算：12回、加減算：3回。その後、正規化するわけですが、掛け算３回、足し算２回、それに、sqrt()という四則演算に比べてズッとコストのかかる計算をしなくてはならないのです。

しかし、球座標で持つようだと、複数回の三角関数の計算に更に時間を要するので結局どちらでも大差なくなってしまいます。

お礼

皆さんありがとうございました.ポイントですが、お二人に同じく20ポイントx10をつけたいぐらいです.特にametsuchiさんの各論の回答はものすごく重要です.今後ともどうぞよろしくお願い致します.

回答No.4

ametsuchi

hogehogeninjaさんの訂正もあったのでお任せしとこうと思いましたが、finetoothcombさんから質問が出てきたので再びシャシャリでてきました。

１）「別方法」：

私は「コストのかかる計算」は嫌いで、外積よりも内積計算を好みます。外積を使って計算時間やコーディングステップ数が減る場合に限って外積を使います。

今の問題の場合、hogehogeninjaさんの言われるように、題意を充たす時、

　n^⊥x^, n^⊥y^　　（１）
だから
n・x^=0、 n・y^=0 （１’）　

です。x^、y^はGiven。「n」は、線上の２パラメータ：α、βに依存しますが、（１’）の条件を課すと、α、βに関する２元一次の連立方程式が得られます。
と言っちゃうと非常にシチメンドクサイことのようですが、紙の上で整理・計算するととても簡単な形になります。

ただ、実質的には、同じ計算量で、外積を使ったからコストがよりかかるという訳ではない（はず）です。ですから、気にしないで下さい。finetoothcombさんにとって分かり易い方を選択してください。

２）「この演算を実行する前段階に、|n_|=0...」

その通りです。２ベクトルの内積による検査を省き、この分母が０に近いかどうかを見ればよいのです。そうすればあとはゼロ割や精度の低下を気にしないでいいです。

ただ、その判定基準、「閾値」だの「トレランス」だのと言いますが、それが実は問題なのです。０割を防ぐだけなら、Cのヘッダファイル（"float.h"）に入っている「FLT_EPSILON」、または「DBL_EPSILON」でもいいのですが、場合によってはもっと大きい数値になる場合もあります。最下層のライブラリとしては、

ａ）引数として外に出すこと。（＝呼び出し側から値を与える）
ｂ）全体として、矛盾を生じないように、「距離」、「角度」など明確な「単位」を持たせることが望ましい。「角度」はコストのかかる計算になりがちなので、正規化したベクトル同士の「内積値」を使う場合がある。

今回の場合、「ベクトル同士の内積値評価は無駄」と言いましたが、それに相当するように、εを手計算で決めておけばいいのです。

上位側、即ちライブラリを使う側を含めてシステム全体で、「ゼロ割りを起こさないよう」気を付けているだけでは、矛盾をきたします。この辺は10分20分で語り尽くせないところなので、省略します。

３）「幾何計算ライブラリ」の目的・利点：

仰せのとおり、「速度向上でなく、再利用性の向上」のためです。処理速度は一般的には遅くなります。

「再利用性の向上」は、広い意味で言われているのだと思いますが、
ａ）開発工数の減少
ｂ）障害発生頻度の減少
ｃ）メンテナンスのし易さ
ｄ）慣れた人がライブラリを開発することにより、品質や、性能（「処理速度」も）が向上
ｅ）後で全く別のシステムに再利用（＝狭義の「再利用性」）

とメリットは限りなくあります。

４）（開発）環境：

Borland C+ 5.5.1 だろうが、VC++だろうが、基本は同じです。AnsiのC++で書けばいいのではないでしょうか？

私は20年以上各種の幾何計算ライブラリの開発に関わってきましたが、基本は、

■処理系に依存しない
■なるべく保守的に書く

ということです。それでも、ここ20年の間になるべく保守的に選んだ言語ですら、

・Fortran４
・Fortran７７
・C（K&R）
・AnsiC
・AnsiC++（VC++だが..）

と変わりました。しかし、上の思想で作っておくと、DOSだろうと、UNIXだろうと、問題なく動きます。一昨年通信会社向けに作ったライブラリもSGIと、PCで問題なく動いています。

しかし去年別の用途で作ったものは心ならずも、

・MFC

を使っています。これだと、SGIやSunなどのEWSに持って行けないです。初めからターゲットマシンが限定されていましたし...。

５）三角関数利用について：

今の持ち方だと、
・直線の通る位置は直交座標
・直線の方向は球座標
という変則的な持ち方ですよね？私の経験では、「球座標」は、

ａ）地球上での位置（緯度・経度）
ｂ）飛行機や自由に３D空間を動くカメラなどの向きを「Euler角」で表現

くらいしか扱ったことがなく、今回の例のように２つの線が「球座標」原点を通る保証がなく一般の捩れの位置にあるというのを大変奇異に感じています。

「球座標」である合理性がなく、基本的に「直交座標」である場合、私の限られた経験では、直線は、
・直線が通る位置（１点）
・直線の方向（単位ベクトル）
で表わすのが普通、というか多かったです。

基本的に「直交座標」でも、カメラ（または「視点」）などを回転させたりするなら、角度や角度を使った三角関数の計算は避けられません。しかし、内部的なデータの持ち方としては、直交座標だけで十分だと私は思います。回転などが加わった場合だけ、座標変換すればいいのです。

三角関数に限定するなら「テーブル云々」は「労多くして益少なし」だと思います。気が付かれているように、角度は任意ですから、補間などしなくてはなりませんから...。

お礼

ありがとうございます.経験のある詳しい人に意見を貰うと、ほんと役立つなぁ、と思いました.

>１）「別方法」：

「"内積"のみ使用して攻める方針は、コスト面で事実上の利がある.例外的な高コスト関数（三角関数等）が並存するなら意味がなくなるかもしれぬが(No.4のご説明)」と理解しました.ものすごく役立ちます.ありがとうございました.なにか、思いつかれるヒントがまた発生されましたら、追加補足していただけましたら幸いです.

>２）「この演算を実行する前段階に、|n_|=0...」
>その通りです。…

"float.h"の、FLT_EPSILON, DBL_EPSILONは初耳でした.あー、すごいなぁ. この***_EPSILONを使わしていただこうと思います.あと(2)の後半のご説明は経験に裏打ちされた深い知識に見えて、簡単には理解しきれなかったのでも少し時間をかけて読もうと思います.コメントありがとうございます.

>３）「幾何計算ライブラリ」の目的・利点：

>処理速度は一般的には遅くなります。
え、そうだったんですか.知らなかったです.意外でした.
> …とメリットは限りなくあります。
こ、これはすごい.なるほど～、かなりメリットですね.

>４）（開発）環境：
>■処理系に依存しない
>■なるべく保守的に書く
なるほど.原則このようにする、といろいろ困ったことも起こらず良いいわけですね.あと、MFCだと移植があまりよろしくないという事実は、とても役立ちます.必要に迫られなければ特に使わないようにしようと思いました.

>５）三角関数利用について：
>…今回の例のように２つの線が「球座標」原点を通る保証がなく一般の捩れの位置にあるというのを大変奇異に感じています。
うーん、これは実に嬉しいご指摘です.なぜなら、「奇異なのか、今私が取り組もうとしている（私が自分に課した）課題は」と解るようになれたからです.現状認識を正しくもてることは、ある意味で、励みになります.特にこのあたりのご指摘はものすごく参考になりました.　いろいろありがとうございます.

回答No.3

hogehogeninja

ametsuchiさんのおっしゃるとおり、
　n^ は x^ × y^
ではなくて、外積を正規化したものです。
このままやると正しい答えになりません。
finetoothcombさんもつられてしまって
n^=a^×b^
と書いてしまっていますが、
　n_=x^×y^
　n^=n_/|n_|
と訂正させてください。

お礼

ありがとうございます.なにしろ知らないことが多いので、なるほど、忘れることがあるから気をつけるべき事なんですね、と、今後に役に立つようで、ありがたい事と感じます.

回答No.2

ametsuchi

hogehogeninjaさんの答えも的確だし、finetoothcombさんもチャンと理解しているようだから、横から口出すのも何ですが...。

１）hogehogeninjaさんの言われる「n^ = x^ × y^ （外積）（２）」の右辺は単位ベクトルとは限らないので、正規化しなくてはなりません。「n^」はそのつもりで書かれたと思うのですが、finetoothcombさん気を付けてください。私は別の方法をよく使いますが、この方法でも計算量は基本的に同じです。

２）hogehogeninjaさんが、「特別な配慮が必要なのは２直線の平行度云々」はとても大事なことです。幾何計算ライブラリやCADシステムを構築する時、一番厄介なのがこの辺です。ただし、「x^ と y^ の内積をとって平行度を調べる」までもなく、（２）で外積計算をして正規化しなくてはならないので、この際自然に検査が必要になります。だから、このやり方なら、「x^ と y^ の内積をとって平行度を調べ」るのは２度手間です。

３）問題自体、「方位角」、「仰角」になっているのでやむをえないのですが、これには三角関数計算が必要です。三角関数はコストの高い計算で、なるべくなら使わない方がよいです。（この場合は仕方ないです！！！）

お礼

ありがとうございます.　有益なコメントに感謝感激です. いろいろな方に方向性の正しさを追認してもらえる、ということにもなり、このジャンルに疎い私としては、ありがたいです.
1-1)外積は正規化が必要なのですね. 気をつけます.
1-2)「別の方法を良く使います(計算量は基本的に同じ）」というのを、お教えいただけたら、とてもありがたいのですが.（どうしても知りたいという訳ではなく、できれば知っておきたいというぐらいなので、もしお時間がありましたらという感じですが.）
2-1)No.3で hogehogeninjaさんが訂正として書いてくださった、n^=n_/|n_|の分母部分ですね.この演算を実行する前段階に、|n_|=0（あるいは、ある小さい値εとして |n_|<ε）かどうかを検査し、もし、それがTRUEなら、平行と見て、処理は中断（というか例外処理）する…という理解をしたのですが、もし間違っていましたら教えてください.
2-2)「幾何計算ライブラリ(やCADシステム）を構築」した経験・知識がなく（やり方が解らない）、Cでべたべたっとプログラムを組むつもりです. 「幾何計算ライブラリ」というのを構築する利点は、速度向上でなく、再利用性の向上というイメージ的理解なのですが、もし間違っていたら教えてくださると幸いです. 使うCの処理系は、DOS上で動くフリーソフトであるところの、Borland C+ 5.5.1 (bcc32.exe)というものを使っているので、それを使う予定です.　これで（仮に余裕があれば経験を積むために）幾何計算ライブラリを構築することはできるのかな.普通どのようにするのでしょうか.
3)そうでした.重要なことを忘れていました.ありがとうございます.コストが高い=計算が遅い、ってことですよね. プログラムを作った後で、結構大量の計算をする予定です.高速化のためのなんか良い作戦がないでしょうか. 値は0≦θ<2π, 0<φ<π/2をとり得るのですが. 一般論としては、高頻度の値については結果の表（テーブル）を用意する作戦がある、と聞いたことがあります.しかし、あらゆる値をとりうる場合には、あまり効果がなさそうです.普通に三角関数入りのプログラムを組むしか仕方ないとは思っています.なにかワザが有りましたら教えていただけましたら幸いです.