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

Question

「空間中の２直線間に、共通垂線を一本引いて、垂線が各直線に対して持つ二つの足の間の距離(---1)と、二つの足の間の中点の位置(---2)」を、求めたいのですが、どのようにするのでしょうか.　

直線は、"平面xy座標上の一点aを通り、かつ、方位角bと仰角cとで定まる方向"という形式で与えられるのです.

最終目標としては、空間中にいろいろと与えられる二直線について、冒頭の「」を求めるプログラムを書きたいのです.

また、冒頭の「」内のことは、捩れの位置で始めて成立するわけですが、交わるのか捩れの位置なのかを、プログラムでもって、簡単に識別する方法ってあるのでしょうか?

このあたりの領域に詳しい方、どうぞ教えてください.

hogehogeninja · Accepted Answer

ベクトルを考えると考えやすいです。

まず、２直線を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^ の内積をとって平行度を調べるなどして）

ametsuchi · Answer

くどいですが、もう一つ補足。

２直線の外積を求めるやり方でも、私流のやり方でもコストは変わらんと言いましたが、実際に計算してみると、私流のやり方の方がコストが安いと分かりました。

直線上のパラメータ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()という四則演算に比べてズッとコストのかかる計算をしなくてはならないのです。

しかし、球座標で持つようだと、複数回の三角関数の計算に更に時間を要するので結局どちらでも大差なくなってしまいます。

ametsuchi · Answer

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角」で表現

くらいしか扱ったことがなく、今回の例のように２つの線が「球座標」原点を通る保証がなく一般の捩れの位置にあるというのを大変奇異に感じています。

「球座標」である合理性がなく、基本的に「直交座標」である場合、私の限られた経験では、直線は、
・直線が通る位置（１点）
・直線の方向（単位ベクトル）
で表わすのが普通、というか多かったです。

基本的に「直交座標」でも、カメラ（または「視点」）などを回転させたりするなら、角度や角度を使った三角関数の計算は避けられません。しかし、内部的なデータの持ち方としては、直交座標だけで十分だと私は思います。回転などが加わった場合だけ、座標変換すればいいのです。

三角関数に限定するなら「テーブル云々」は「労多くして益少なし」だと思います。気が付かれているように、角度は任意ですから、補間などしなくてはなりませんから...。

hogehogeninja · Answer

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

ametsuchi · Answer

hogehogeninjaさんの答えも的確だし、finetoothcombさんもチャンと理解しているようだから、横から口出すのも何ですが...。

１）hogehogeninjaさんの言われる「n^ = x^ × y^ （外積）（２）」の右辺は単位ベクトルとは限らないので、正規化しなくてはなりません。「n^」はそのつもりで書かれたと思うのですが、finetoothcombさん気を付けてください。私は別の方法をよく使いますが、この方法でも計算量は基本的に同じです。

２）hogehogeninjaさんが、「特別な配慮が必要なのは２直線の平行度云々」はとても大事なことです。幾何計算ライブラリやCADシステムを構築する時、一番厄介なのがこの辺です。ただし、「x^ と y^ の内積をとって平行度を調べる」までもなく、（２）で外積計算をして正規化しなくてはならないので、この際自然に検査が必要になります。だから、このやり方なら、「x^ と y^ の内積をとって平行度を調べ」るのは２度手間です。

３）問題自体、「方位角」、「仰角」になっているのでやむをえないのですが、これには三角関数計算が必要です。三角関数はコストの高い計算で、なるべくなら使わない方がよいです。（この場合は仕方ないです！！！）

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

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

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

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

ametsuchiさんのおっしゃるとおり、

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

似たような質問が見つかりました

関連するカテゴリからQ&Aを探す

デイリーランキングこのカテゴリの人気デイリーQ&Aランキング

マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング