3次元空間上の点A(x1, y1, z1)と点B(x2, y2 z2)を結んで出来る線分の中点を知りたいのですが、完全な文系出身であまり数学に詳しくないため、公式の見方がよくわかりません。 Wikipediaの中点のページにあるn次元ユークリッド空間上の中点の公式がそれのようですが、「n 次元ユークリッド空間上の2点 A, B を直交座標系であらわし、それぞれを (a0, ..., an－1), (b0, ..., bn－1) とすると」の時点ですでに理解できないので、単純な公式で教えて下さると助かります。

中点座標を(x3,y3,z3)とすると x3=(x1+x2)/2 y3=(y1+y2)/2 z3=(z1+z2)/2

3次元空間上の2点を結ぶ線分の中点を知りたい -3次元空間上の点A(x1, y- 数学

Q3次元座標2点からの直線式の求め方

お世話になります。

3次元座標2点からの直線式(ax+by+cz=0)の求め方を教えて下さい。

2次元座標であれば、1つの傾きから算出できるのですが、3次元座標になると、X-Y平面、Y-Z平面での傾きの使い方がこんがらかってしまいます。
基本的な質問で申し訳ありませんが、よろしくお願い致します。

座標1 = (x1,y1,z1)
座標2 = (x2,y2,z2)

以上

Aベストアンサー

> 直線式(ax+by+cz=0)の求め方を教えて下さい。
3次元座標では(ax+by+cz=0)は原点を通る平面になり、直線の式ではありません。ax+by+cz=dは平面の一般式です。

2点を通る直線の式には公式があります。
以下のように簡単に導けます。
点(x1,y1,z1)を通り方向ベクトル(x2-x1,y2-y1,z2-z1)の直線ですから
媒介変数形式で
(x,y,z)=(x1,y1,z1)+t(x2-x1,y2-y1,z2-z1)
と成ります。
これを変形してすれば
(x-x1)/(x2-x1)=(y-y1)/(y2-y1)=(z-z1)/(z2-z1)
と3次元座標の直線の式となります。

他の回答も見る

Q３次元空間上の２つの座標から角度を求めたい

３次元空間上の２つの座標、座標A と座標B から
A -> B の角度を求めるにはどうしたらよいでしょうか？

Pythonでプログラムを組んでいるのですが調べてもイマイチ３次元での方法がわからなかったため質問しました。

多分座標Aにオブジェクトがあると仮定した時、
そのオブジェクトを座標Bに向けるというのと同じようなことをすればいいのではと考えて
調べているのですがわかりません。

２次元ではatanを利用すれば出来るようなのですが
３次元の方法がわからなかったため質問しました。

座標A：(xA, yA, zA)
座標B：(xB, yB, zB)

求めたい結果：Rotation(xR,yR,zR)

空間的には

+X 左
+Y 上
+Z 正面

になっています。

Pythonの方法でなくても言語などは問わないので分かる人教えて下さい。
よろしくお願いします。

Aベストアンサー

あなたが期待しているのは、どんな結果なのですか?
「これで近い向きにはなっていますがやはりずれています」とは、どんな方法で確認したものなのでしょうか?

A(5,5,5) から B(10,10,10)への角度は
B-A=(10-5,10-5,10-5) = (5,5,5) のベクトルと並行な角度になります。
角度をこのままX,Y,Z座標で表す方法もあります。
大きさで割った単位ベクトルを角度として使用することもあります。

平面なら角度は1つで済みますが、立体では2つ必要です。
X-Y平面上でどの向きかを示す「方位角」と上下にどのくらいかの「仰角」の組合せはよく使われます。

(5,5,5)ベクトルの方位角は45度、仰角は35.26...度で、合ってます。

何に使ってるかわかりませんが、Vector3D[0]が方位角、Vector3D[2]に仰角、って使い方が間違っているということは無いですか?

他の回答も見る

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

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

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

Aベストアンサー

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

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

３）で、以下、２）のチェック済みである...続きを読む

他の回答も見る

Q円の中心座標の問題の解き方を教えてください。

円の中心座標の問題の解き方を教えてください。

問題
3点A(1,-2,1) B(3,1,7) C(2,0,6)を通る円の中心座標を求めよ

解く方針は以下のようにしました。
中心座標をO(X,Y,Z)とおき、AB、BCとCAの中点をそれぞれ点D,E,Fとして
OD,OE,OFはそれぞれAB、BC、CAと直交することから内積を利用して中心Oを求めようとしましたが、
解けませんでした。

よろしくお願いします。

Aベストアンサー

3点が１平面上にあって、これらを通る円というのは気がつきませんでした。
3点A(1,-2,1) B(3,1,7) C(2,0,6)を通る円はこの平面では円でもｘｙ平面に投影すると楕円になるところがにくいところです。

（解答)
円の中心座標をP(x,y,z)、半径をrとすると
(x-1)^2+(y+2)^2+(z-1)^2=r^2 (1)
(x-3)^2+(y-1)^2+(z-7)^2=r^2　　(2)
(x-2)^2+y^2+(z-6)^2=r^2　　　　(3)

(2)－(3)より
x+y+z=19/2 (4)
(1)－(3)より
x+2y+5z=17　　(5)

点P(x,y,z)、A(1,-2,1)、 B(3,1,7)、 C(2,0,6)
が1平面状にあることから

行列式
|...続きを読む

他の回答も見る

Q二点の座標から角度を求めるには？

2点の座標A,Bの角度を求めたいのですが,たとえばA点(0,0)とB点(4,3)を結ぶラインは、底辺Bxと高さByを元に三角関数？から30度と求められますが、B点がマイナス座標が絡んできた場合などの90度から359度までをどう求めていいか悩んでいます。また、A点も(0,0)に限定されるわけではないので、ますます混乱しています。どう考えればよいのか教えていただきたいのですが
（水平はX軸プラス方向が0度です）

Aベストアンサー

>2点の座標A,Bの角度を求めたい～・・・・

このままなら答えは０ですけど？

ｘｙ座標で、ｘ軸のプラス方向を０度とし、
２点の座標Ａ、Ｂにより形成される線ＡＢとｘ軸との角度
ってことですね。

>たとえばA点(0,0)とB点(4,3)を結ぶラインは、底辺Bxと高さByを
>元に三角関数？から30度と求められますが、

sen-senさんの書かれたとおり、これは間違いです。
この場合、Ｂからｘ軸へのばした垂線とｘ軸との交点をＣとすると、
三角形ＡＢＣができ、そのときの求めたい角度をθとすると、
tanθ=3/4...続きを読む

他の回答も見る

QLNK2019: 未解決の外部シンボルのエラーが出る

Microsoft Visual Studio 2008
Version 9.0.21022.8 RTM
Microsoft .NET Framework
Version 3.5 SP1
----------------------------------------------------------------
新しいプリジェクト→Win32 コンソール　アプリケーション（ソリューションのディレクトリを作成　チェック外す）→Windows　アプリケーション（空のプロジェクト　チェック外す）
----------------------------------------------------------------
　プログラム

　mymain.cpp
#include "myhelper.h"
#include "mymain.h"

//自キャラのデータ
Point2D g_jikipos = {40, 400};//自キャラの座標

//画像ハンドル
int g_jikiimage[11];

//色々なファイルの読み込み
int LoadFiles(){
//画像ファイル読み込み
if(LoadDivGraph("media\\player01.bmp",
11,11,1,64,64,g_jikiimage) == -1) return -1;

return 1;
}

　mymain.h
//他から呼び出させるMyMainの関数
void MyMain();
int LoadFiles();

　myhelper.h（サンプルなので打ちミスはない）
#include "DxLib.h"
#include <limits.h>
#include <math.h>

//構造体宣言
//座標またはベクトルを記録する構造体
struct Vector{
float x,y;
};
typedef Vector Point2D;
//線を記録する構造体
struct Line2D{
Point2D startpos, endpos;
float katamuki;//傾きをラジアン値で記録
Vector speed;//移動している場合は速度をセット
};
//球体を記録する構造体
struct Ball2D{
Point2D position;
float hankei;//半径
};
//四角形を記録する構造体
struct Rect2D{
Point2D lefttop;
Point2D rightbottom;
float width;
float height;
};

//ライブラリ関数
Point2D PosInView(Point2D in);
int XInView(float inx);
int YInView(float iny);
void ScrollToLeft(float jikiposx);
void ScrollToRight(float jikiposx);
void ScrollToUp(float jikiposy);
void ScrollToDown(float jikiposy);
void DrawLineInView(float x1, float y1, float x2, float y2, int Color, int Thickness);
void DrawCircleInView(float x, float y, float r, int Color, int FillFlag);
void DrawAnimation(float x, float y, double ExtRate, double Angle,int TurnFlag,
int *imgarray, int allframe, float fps);
//ベクトル関数
Vector CreateVector(Vector in, float veclen);
Vector AddVector(Vector v1, Vector v2);
Vector SubVector(Vector v1, Vector v2);
Vector AddVectorInFrameTime(Vector pos, Vector speed);
Vector AddVectorInFrameTime2(Vector pos, Vector speed, Vector accel);
Vector Normalize(Vector in);
Vector RotateVector(Vector in, float radian);
float VectorLengthSquare(Vector in);
float DotProduct(Vector v1, Vector v2);
float CrossProduct(Vector v1, Vector v2);
void SetLine2DKatamuki(Line2D *in);
void DrawLine2D(Line2D in, int Color, int Thickness);
void DrawBall2D(Ball2D in, int Color, int Fill);
//当たり判定関数
bool HitTestLineAndBall(Line2D linein, Ball2D ballin);
bool IsPointAtLineFace(Line2D linein, Point2D ptin);
bool HitTestLineAndLine(Line2D line1, Line2D line2);
bool HitTestBallAndBall(Ball2D a, Ball2D b);
bool HitTestPointAndBox(Rect2D rect, Point2D pt);
//タイマー関数
void SetSimpleTimer(int idx, int time);
int GetPassedTime(int idx);

//グローバル変数
extern float g_frametime;
extern Rect2D g_framerect;//画面領域（当たり判定）
extern Point2D g_current_field_pos;//現在の左上座標
extern Rect2D g_stagesize;//ステージサイズ

//定数宣言
const float ZEROVALUE = 1e-10f;
const float PIE = 3.1415926f;
const int SCROLL_LIMIT = 200;
----------------------------------------------------------------
　エラー内容
1>myhelper.obj : error LNK2019: 未解決の外部シンボル "void __cdecl MyMain(void)" (?MyMain@@YAXXZ) が関数 _WinMain@16 で参照されました
1>C:\Documents and Settings\Owner\My Documents\Visual Studio 2008\Projects\my\Debug\my.exe : fatal error LNK1120: 外部参照 1 が未解決です
1>my - エラー 2、警告 0
ビルド: 0 正常終了、1 失敗、0 更新不要、0 スキップ
----------------------------------------------------------------
画像を貼り付けときます
（見えにくい場合→http://www.dotup.org/uploda/www.dotup.org154142.jpg.html）
初心者なのでわかりやすくお願いします

Microsoft Visual Studio 2008
Version 9.0.21022.8 RTM
Microsoft .NET Framework
Version 3.5 SP1
----------------------------------------------------------------
新しいプリジェクト→Win32 コンソール　アプリケーション（ソリューションのディレクトリを作成　チェック外す）→Windows　アプリケーション（空のプロジェクト　チェック外す）
----------------------------------------------------------------
　プログラム

　mymain.cpp
#include "myhelper.h"
#include "mymain.h"

//自...続きを読む

Aベストアンサー

ファイル構成から推測するに
mymain.cpp というファイルに
void MyMain(void) {
// ここに処理を書く
}
という関数が必要なようです。

他の回答も見る

Q三点の座標から中心点の求め方

三点の座標から中心点を求める公式は以下のHPからわかったのですが、なぜこのような式が成り立つのかがわかりません。

http://www.geocities.co.jp/SiliconValley-Oakland/1253/3point.html

分かる方がおられましたらお教え頂けないでしょうか。
参考URLがありましたらそちらでも構いません。
宜しくお願いします。

Aベストアンサー

表記されているサイトの考え方では三点を周上に含む円の中心を求めています。
円の定義というのはある１点から距離の等しい点の集まりです。
そこから考えれば円の中心点というのは、円周上の全ての点から距離の等しい点ということになります。
なので、中心点を求めるということは三点から距離の等しい点を求めればよいことになります。
三点を取れば一点に決まるためそれ以上は取らなくてもかまわないからです。

二点を結ぶ線分の垂直二等分線が二点間から距離の等しい点の集まりですから、サイトの表記方法に倣うとP1とP2の垂直二等分線とP2とP3の垂直二等分線の交点はP1,P2,P3の三点からの距離が等しい、すなわち三点を円周上に含む円の中心ということになります。

※質問者の方のレベルがどの程度かわからないのでこのような簡易的な回答にしました

他の回答も見る

Q平面方程式の傾きについて

ax+by+c=zの平面方程式の
傾きというものは存在するものなのですか？

もし存在するのであればどなたか
教えていただけないですか？

Aベストアンサー

こんばんは。

ｙ　＝　ａｘ　＋　ｂ
これの傾きはご存知かと思います。
しかし、「傾き」と言うからには、「何に対しての傾き」かという定義が必要ですよね。
上記の式の傾きは、Ｘ軸に対する傾きです。
Ｘ軸に対する傾きを角度θを用いて表せば、
Δｙ　＝　Δｘ・tanθ
です。（tanθ　＝　ａ）
そして、ｘに対するｙの傾きも考えることもできます。
それは、当然ながら、１／ａ　です。

さて、
同様に、平面の方程式を考えるときにも、何に対しての傾きを求めるかを決めなければいけません。
それは３通りあります。
・Ｘ－Ｙ平面　（平面ｚ＝０　に同じ）
・Ｙ－Ｚ平面　（平面ｘ＝０　に同じ）
・Ｚ－Ｘ平面　（平面ｙ＝０　に同じ）

平面と平面との傾きを求めるということは、それらの法線同士の傾きを求めることと同じです。

というわけで、平面の法線の方程式の求め方を学んでください。

他の回答も見る

Q3次元ベクトルをある軸ベクトルで回転させたい

3次元ベクトルの求め方を教えてください。

下記図のように始点を軸ベクトルでθ（度）だけ回転したときの？の位置を求めたいのです。
これはどのような計算方法になるのでしょうか？なかなか思いつかなくて非常に悩んでいます。
アドバイスや回答をいただけたら助かります。よろしくお願いします。

Aベストアンサー

先ず、中心点(Sx,Sy,Sz)が原点にくるよう全体を平行移動させます。
（一番最後に元に戻します）
始点(Px,Py,Pz)は、(Px-Sx,Py-Sy,Pz-Sz)に移ります。この座標を(Px',Py',Pz')とします。

次に、回転軸ベクトル（Ax Ay Az)を回転させ、x軸に合致させます。それには二回の
回転変換が必要です。
最初に、ベクトル（Ax Ay Az)と、x軸方向単位ベクトル(1 0 0）のなす平面の法線ベクトルが
z軸に合うよう、x軸を回転させます（その角度をφとします）。
すると、回転軸ベクトルはx-y平面上に乗るので、それがx軸に合うよう、z軸を回転させます
（その角度をψとします）。

ベクトル（Ax Ay Az)と、x軸方向単位ベクトル(1 0 0）のなす平面の法線ベクトルは、(0 Az -Ay)。
x軸周りにφ回転させると、このベクトルは、
「1　　0　　　　0　　　「 0　　=「　　　　　　0
0　cosφ　-sinφ 　 Az　　　Az・cosφ+Ay・sinφ
0　sinφ　 cosφ」　-Ay」　　Az・sinφ-Ay・cosφ」
で、z軸ベクトルに合うので
「　　　　　　0　　　　　 =「0
Az・cosφ+Ay・sinφ　　0　
Az・sinφ-Ay・cosφ」　 1」
これから、cosφ=-Ay/(Ay^2+Az^2)、sinφ=Az/(Ay^2+Az^2)
∴　φ=Arctan(-Az/Ay)

回転軸ベクトル（Ax Ay Az)は、
「1　　0　　　　0　　　「Ax　=「　　　　　　Ax　　　　　　=「　　　　　　　Ax　　　　　　　　　　　　　　　　　　　=「Ax　
0　cosφ　-sinφ 　 Ay　　　Ay・cosφ-Az・sinφ　　 Ay・{-Ay/(Ay^2+Az^2)}-Az・{Az/(Ay^2+Az^2)}　　 -1
0　sinφ　 cosφ」　 Az」　　 Ay・sinφ+Az・cosφ」　 Ay・{Az/(Ay^2+Az^2)}+Az・{-Ay/(Ay^2+Az^2)}」　　0」
に変換され、x-y平面上に乗ります。これを(Ax' Ay' Az') とします。
つまり、(Ax' Ay' Az')=(Ax -1 0)

始点(Px',Py',Pz')もこの変換を受けるのですが、変換を全部纏めて後、一括変換させます。

今度は、x-y平面上に乗った回転軸ベクトル（Ax' Ay' Az')を、z軸の周りにψ回転させます。
「cosψ　-sinψ　0　「Ax' 　=「Ax'・cosψ-Ay'・sinψ　=「Ax・cosψ+sinψ
sinψ　 cosψ　0 　 Ay'　　 Ax'・sinψ+Ay'・cosψ　　 Ax・sinψ-cosψ
　 0　　　　0　　 1」　 Az'」　　　　　　　Az'　　　　　　」　　　　　0　　　　　　」
これが、x軸ベクトルに合うので、
Ax・cosψ+sinψ=1
Ax・sinψ-cosψ=0
これから、cosψ=Ax/(Ax^2+1)、sinψ=1/(Ax^2+1)
∴　ψ=Arctan(1/Ax)

以上の回転の変換の積は、
「cosψ　-sinψ　0　「1　　0　　　　0　　 =「cosψ　-sinψ・cosφ　 sinψ・sinφ
sinψ　 cosψ　0 　 0　cosφ　-sinφ 　 sinψ　　cosψ・cosφ　-cosψ・sinφ
　 0　　　　0　　 1」　 0　sinφ　 cosφ」　　 0　　　　　sinφ　　　　　cosφ　　　」

この変換を始点(Px',Py',Pz')に施します。
「cosψ　-sinψ・cosφ　 sinψ・sinφ　　「Px' = 「Px'・cosψ-Py'・sinψ・cosφ+Pz'・sinψ・sinφ
sinψ　　cosψ・cosφ　-cosψ・sinφ　 Py'　　 Px'・sinψ+Py'・cosψ・cosφ-Pz'・cosψ・sinφ
　 0　　　　　sinφ　　　　　cosφ　　　」　Pz'」　　Py'・sinφ+Pz'・cosφ　　　　　　　　　　　　　　」　

この点を(Px”,Py”,Pz”)とします。

さて、ここでx軸に合った回転軸ベクトル(1 0 0)周りに(Px”,Py”,Pz”)を角度θ、回転させます。
「1　　0　　　　0　　　「Px”　=「　　　　　Px”　　　
0　cosθ　-sinθ 　 Py”　　Py”・cosθ-Pz”・sinθ　
0　sinθ　 cosθ」　 Pz”」　 Py”・sinθ+Pz”・cosθ」

これを(P_x, P_y, P_z）とします。

今度は、回転させた回転軸を元に戻す変換です。
回転の変換の逆行列は、行列各要素の余因子の行と列を入れ替えたものを行列式で割ったもので、
行列式は、(cosψ）^2+(sinψ)^2=1 なので、逆行列は
「　cosψ　　　　　　sinψ　　　　　　　0　　
-sinψ・cosφ　　cosψ・cosφ　　 sinφ
sinψ・sinφ　　 -cosψ・sinφ　　cosφ」

これを、(P_x, P_y, P_z）に施します。
「　cosψ　　　　　　sinψ　　　　　　　0　　「P_x　=「P_x・cosψ+P_y・sinψ
-sinψ・cosφ　　cosψ・cosφ　　 sinφ　　P_y　　 -P_x・sinψ・cosφ+P_y・cosψ・cosφ+P_z・sinφ
sinψ・sinφ　　 -cosψ・sinφ　　cosφ」 P_z」　 P_x・sinψ・sinφ-P_y・cosψ・sinφ+P_z・cosφ」

結局、θ回転後のP点の座標は、
x座標： P_x・cosψ+P_y・sinψ
y座標： -P_x・sinψ・cosφ+P_y・cosψ・cosφ+P_z・sinφ
z座標： P_x・sinψ・sinφ-P_y・cosψ・sinφ+P_z・cosφ
となります。

ここで、置き換えた変数を順次、元に戻します。
P_x、P_y、P_z を Px”、Py”、Pz” に、
Px”、Py”、Pz” を Px’、Py’、Pz’ に、
最後に、平行移動を戻して　Px’、Py’、Pz’ を Px、Py、Pz に直します。

先ず、中心点(Sx,Sy,Sz)が原点にくるよう全体を平行移動させます。
（一番最後に元に戻します）
始点(Px,Py,Pz)は、(Px-Sx,Py-Sy,Pz-Sz)に移ります。この座標を(Px',Py',Pz')とします。

次に、回転軸ベクトル（Ax Ay Az)を回転させ、x軸に合致させます。それには二回の
回転変換が必要です。
最初に、ベクトル（Ax Ay Az)と、x軸方向単位ベクトル(1 0 0）のなす平面の法線ベクトルが
z軸に合うよう、x軸を回転させます（その角度をφとします）。
すると、回転軸ベクトルはx-y平面上に乗るので、それがx軸...続きを読む

他の回答も見る

Q３次元座標上の２直線の交点判定について

座標A(x1,y1,z1)から座標B(x2,y2,z2)への線分ABと
座標C(x3,y3,z3)から座標B(x4,y4,z4)への線分CDがあり、
線分ABと線分CDが交点を持つかどうかのプログラムを作りたいです。

C言語かVBかFortranで記述され、DirectXやOpenGLのライブラリを使わない方法の
サンプルソースの載っているページを教えていただけませんか？
また、ご迷惑でなければソースコードを記述していただけると助かります。

Aベストアンサー

2線分の最短距離を求めてそれが0になる場合，交わっていると判断するのが良いと思います．
piyoは2線分の最短距離を求めるサブルーチンです．
2線分が平行なときは別の処理が必要ですね．

subroutine piyo(r1,r2,r3,r4,d)
implicit none
real(8),intent(in)::r1(3),r2(3),r3(3),r4(3)
real(8),intent(out)::d
real(8) r12(3),r34(3),b(2),a(2,2),s,t,rr(3),deta
r12(:)=r2(:)-r1(:)
r34(:)=r4(:)-r3(:)

!!! 2つの直線の距離の2乗をdとする
!!! d=|((r1+s*r12)-(r3+t*r34))|^2
!!! dの最小値を求める
!!! dをsについて偏微分
!!! r12・((r1+s*r12)-(r3+t*r34))=0
!!! dをtについて偏微分
!!! -r34・((r1+s*r12)-(r3+t*r34))=0
!!! sとtについてまとめると
!!! r12・r12 s - r12・r34 t = -r12・r1 +r12・r3
!!! -r34・r12 s + r34・r34 t = r34・r1 -r34・r3
a(1,1)=dot_product(r12,r12)
a(2,1)=-dot_product(r12,r34)
a(1,2)=a(2,1)
a(2,2)=dot_product(r34,r34)
b(1)=-dot_product(r12,r1)+dot_product(r12,r3)
b(2)=dot_product(r34,r1)-dot_product(r34,r3)
deta=a(1,1)*a(2,2)-a(1,2)*a(2,1)
!!! sとtについて解く
if(deta /=0.d0) then
s=(a(2,2)*b(1)-a(1,2)*b(2))/deta
t=(-a(2,1)*b(1)+a(1,1)*b(2))/deta
!!! elseif
!!! deta=0の場合は別途処理する必要あり．2つの線分が平行な場合deta=0となる．
end if

!!! この時点で0<s<1, 0<t<1でなければ，交わっていないことが分かる．
!!! if(s<0.d0.or.1.d0<s.or.t<0.d0.or.1.d0<t) then
!!! write(*,*) "交わってない"
!!! stop

if(s<0.d0) then
s=0.d0
elseif(1.d0<s)then
s=1.d0
end if
if(t<0.d0) then
t=0.d0
elseif(1.d0<t)then
t=1.d0
end if

!!! 最短距離の計算
rr(:)=r1(:)+s*r12(:)-(r3(:)+t*r34(:))
d=sqrt(dot_product(rr,rr))
end subroutine piyo

program hoge
implicit none
real(8) r1(3),r2(3),r3(3),r4(3),d
r1(:)=(/0.d0,0.d0,0.d0/)
r2(:)=(/1.d0,1.d0,1.d0/)
r3(:)=(/1.d0,0.d0,0.d0/)
r4(:)=(/1.d0,1.d0,0.d0/)
call piyo(r1,r2,r3,r4,d)
if(d<1.d-15) then
write(*,*) "交わる"
endif
end program

2線分の最短距離を求めてそれが0になる場合，交わっていると判断するのが良いと思います．
piyoは2線分の最短距離を求めるサブルーチンです．
2線分が平行なときは別の処理が必要ですね．

subroutine piyo(r1,r2,r3,r4,d)
implicit none
real(8),intent(in)::r1(3),r2(3),r3(3),r4(3)
real(8),intent(out)::d
real(8) r12(3),r34(3),b(2),a(2,2),s,t,rr(3),deta
r12(:)=r2(:)-r1(:)
r34(:)=r4(:)-r3(:)

!!! 2つの直線の距離の2乗をdとする
!!! d=|((r1+s*r12)-(r3+t*r34))|^2
!!! dの最小値を求め...続きを読む