No.4ベストアンサー
- 回答日時:
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
この回答への補足
質問者のGyustabです。
以下3点質問させてください。
(1)!!! 2つの直線の距離の2乗をdとする
!!! d=|((r1+s*r12)-(r3+t*r34))|^2
!!! dの最小値を求める
の部分について
線分r1r2上の点を(r1+s*r12)
線分r3r4上の点を(r3+t*r34) とすると
2点間の距離は
sqrt( ( (r3+t*r34)-(r1+s*r12) )^2 ) ⇒ (r3+t*r34)-(r1+s*r12) となり
距離の2乗は ((r3+t*r34)-(r1+s*r12))^2 となるという理解で良いでしょうか?
(2)距離の2乗をdとしているのは、
s,tについて微分出来る形にするためわざと2乗していると考えれば良いですか?
(3) 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)
上記数式がわかりません、
一連の流れは内積を使って何を計算しているのでしょうか?
わざわざFortranのプログラムを記述いただきありがとうございます。
本日、線分が半径を持つことが判明したため、
線分が衝突しないためには、
線分間で 半径*2 より大きい離隔距離が必要になりました
そのため nineexitさんのコードを解析し、理解したく思います。
まず、
!!! 2つの直線の距離の2乗をdとする
!!! d=|((r1+s*r12)-(r3+t*r34))|^2
の部分がわからず、なぜこのように定義できるのか現在ネットで調査しています。
dot_product はFortranの関数でベクトル間の内積を求めるということがわかりました、
以下のページを参考にエクセルVBAのコードに落とせるよう調べています。
http://homepage1.nifty.com/gfk/vector-naiseki.htm
まだまだ時間がかかりますが少しづつコードを解析していきます。
No.6
- 回答日時:
> つまり、2直線間の距離の2乗であるdは、
> 2直線間の距離同士の内積と考えるということでしょうか?
2つの直線上にある"2つの点"の距離の自乗がdです.
「2直線間の距離同士の内積」ではありません.
内積は距離同士でとるものではなくベクトル同士でとるものです.
同じベクトル同士の内積を取ると長さの自乗が得られるので,内積を取っています.
>2つの直線上にある"2つの点"の距離の自乗がdです.
>「2直線間の距離同士の内積」ではありません。
>内積は距離同士でとるものではなくベクトル同士でとるものです。
>同じベクトル同士の内積を取ると長さの自乗が得られるので,内積を取っています.
ご指摘ありがとうございます、数学素人の回答ですいません、
ベクトルと距離がごちゃまぜになっていました、
ご指摘いただいた点および今まで教えていただいた点を
含めて一度ソースコードを納得できるところまでコーディングしてみます。
No.5
- 回答日時:
回答に対する質問への回答です。
(1)
r1, r2, r3, r4, r12, r34はベクトルです。
ベクトルa=(a1,a2,a3)の長さの自乗はa同士の内積(a・a)になります。
サブルーチンの終わりに
rr(:)=r1(:)+s*r12(:)-(r3(:)+t*r34(:))
d=sqrt(dot_product(rr,rr))
とありますが、これが長さの自乗を計算しているところです。
(2)
微分しやすいように自乗しています。
2直線上の点の距離の最小値を求めるためには、2直線上の点の距離の自乗の最小値を求めればよいので、このようなことをしています。
最後に最短距離を算出するときは、きちんとルートを取っています。
(3)
!!! r12・r12 s - r12・r34 t = -r12・r1 +r12・r3
!!! -r34・r12 s + r34・r34 t = r34・r1 -r34・r3
を行列A、ベクトルbを用いて表記しました。
Ax=bの形でxはsとtからなるベクトルです。x=(s,t)
detaは行列A行列式です。
ただの連立一次方程式なので、行列やベクトルを持ち出す必要はありません。
返答ありがとうございます。
>ベクトルa=(a1,a2,a3)の長さの自乗はa同士の内積(a・a)になります。
つまり、2直線間の距離の2乗であるdは、
2直線間の距離同士の内積と考えるということでしょうか?
>d=sqrt(dot_product(rr,rr))
>とありますが、これが長さの自乗を計算しているところです。
ベクトルa=(a1,a2,a3)の長さの自乗はa同士の内積(a・a)ということなので、
dot_product(rr,rr)はベクトルrrの長さの2乗と同じ意味ですね、
ということはベクトルrr同士の内積の平方根を求めるとベクトルrrの長さが出てくるという事ですね
>微分しやすいように自乗しています。
>2直線上の点の距離の最小値を求めるためには、
>2直線上の点の距離の自乗の最小値を求めればよいので、このようなことをしています。
>最後に最短距離を算出するときは、きちんとルートを取っています。
とてもわかりやすかったです、
紙上に数式を記述し、展開してs,tに対して微分したところ、
sとtについてまとめた数式まで辿り着けました。
No.3
- 回答日時:
#1です。
>下記4つの数式の条件を満たすのは4点の座標の値が整数の場合のみでしょうか?
>座標の値が単精度浮動小数点や倍精度浮動小数点の場合は有効桁数の範囲内であれば成立しますか?
数学的には、実数の範囲で成立します。
コンピュータで単精度、倍精度で計算する場合は当然有効桁数の問題がでてきます。
特に、減算で有効桁数が桁落ちする可能性がでてくるので、十分考慮する必要があります。
前回の回答で、最初に「コードを書くのは面倒」と書いたのはそういう意味です。悪しからず。
ご返答いただきありがとうございます。
>前回の回答で、最初に「コードを書くのは面倒」と書いたのはそういう意味です。
>悪しからず。
いえいえ、疑問に対する答えがわかっただけで助かりました。
明日、会社にて数式の許容誤差や
計算上の有効桁数の桁落ちをどのようにすれば良いかを熟考してみます。
No.2
- 回答日時:
学校の宿題か何かの回答をここに丸投げ質問してないですか。
そういうのは先般の大学入試問題にでも批判を浴びたはず。もしそうなら、授業で先生に教えてもらえるのではないですか。それをここで教えるのは授業妨害とも言える。
Googleででも「空間 2直線 交点座標 」などで照会すれば記事はある。
ベクトルや行列の知識が必要な記事が多いが。そういうのは理解できるのか。
http://www.teu.ac.jp/aqua/GS/text/PDF/3DMath.pdf
>学校の宿題か何かの回答をここに丸投げ質問してないですか。そういうのは先般の大学入試問題にでも>批判を浴びたはず。
>もしそうなら、授業で先生に教えてもらえるのではないですか。それをここで教えるのは授業妨害とも>言える。
いいえ違います社会人です、仕事上必要であったため質問しました、
教えてくれる先生などいませんし、授業を妨害している暇はありません。
>Googleででも「空間 2直線 交点座標 」などで照会すれば記事はある。
「線分 交点 3次元」、「直線 衝突判定 3次元」、「直線 交点 サンプルプログラム」等で
検索しましたが、概念はたくさんあっても私の解読できるC言語・VisualBasic・Fortran・javaのソースが発見できなかったためここに質問しました。
客先にはエクセルで提供するのでDirectXやOpenGLのようなライブラリを要するものは使用したくないのです。
>ベクトルや行列の知識が必要な記事が多いが。そういうのは理解できるのか。
数学者でも教師でも数値解析エンジニアでもないのでベクトルや行列について1~100まで
理解する気はありません、客先の仕様を満たすために必要な概念と考え方、
ソースコードへ落とす事が必要なだけです。
以上、なにか返答があればどうぞ
No.1
- 回答日時:
コードを書くのは面倒なので考え方だけ。
x12=x2-x1、y12=y2-y1、z12=z2-z1、
x34=x4-x3、y34=y4-y3、z34=z4-z3とすると、
線分AB上の点は(x1+s*x12, y1+s*y12, z1+s*z12) (0≦s≦1)
線分CD上の点は(x3+t*x34, y3+t*y34, z3+t*z34) (0≦t≦1)
と表現できます。
線分ABと線分CDが交点を持つとすると、
(1) x1+s*x12=x3+t*x34
(2) y1+s*y12=y3+t*y34
(3) z1+s*z12=z3+t*z34
が成り立ちます。
(1),(2)式からs,tを求めると、
s=((x3-x1)*y34-x34*(y3-y1))/(x12*y34-x34*y12)
t=((x1-x3)*y12-x12*(y1-y3))/(x34*y12-x12*y34)
このS,tが0≦s≦1 , 0≦t≦1 の条件を満たし、かつ、
このs,tを(3)式に代入して等号が成立するときに交点を持ちます。
あとは、これをコードにすればいいだけです。
補足
x12*y34=x34*y12 の場合、2直線は、同じか、並行か、ねじれ位置にあるので交差していないことになります。
この回答への補足
質問者のGyustabです。
コーディング中疑問に思ったため再度質問させてください。
下記4つの数式の条件を満たすのは4点の座標の値が整数の場合のみでしょうか?
座標の値が単精度浮動小数点や倍精度浮動小数点の場合は有効桁数の範囲内であれば成立しますか?
成立しないのであれば (x12*y34)/(x34*y12)≒1.00000 のように限りなく1に近づく誤差を
計算したり、s=・・・・やt=・・・の式の乗算・除算の順番を考え有効桁数の
精度を考慮した計算をする必要があるのではないかと思っています。
数式1:s=((x3-x1)*y34-x34*(y3-y1))/(x12*y34-x34*y12)
数式2:t=((x1-x3)*y12-x12*(y1-y3))/(x34*y12-x12*y34)
数式3:z1+s*z12=z3+t*z34
数式4:x12*y34=x34*y12
返信ありがとうございます、以下の部分が大変助かりました。
s=((x3-x1)*y34-x34*(y3-y1))/(x12*y34-x34*y12)
t=((x1-x3)*y12-x12*(y1-y3))/(x34*y12-x12*y34)
さまざまなホームページで紹介されている概念の通りにやれば紙の上では解けるのですが、
どうやって2線分の交点をソースコードに落とせば良いかで四苦八苦しておりました。
今よりコーディングしてみます。
補足もありがとうございます、この部分も条件に組み込みます。
お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!
関連するカテゴリからQ&Aを探す
おすすめ情報
- ・漫画をレンタルでお得に読める!
- ・人生のプチ美学を教えてください!!
- ・10秒目をつむったら…
- ・あなたの習慣について教えてください!!
- ・牛、豚、鶏、どれか一つ食べられなくなるとしたら?
- ・【大喜利】【投稿~9/18】 おとぎ話『桃太郎』の知られざるエピソード
- ・街中で見かけて「グッときた人」の思い出
- ・「一気に最後まで読んだ」本、教えて下さい!
- ・幼稚園時代「何組」でしたか?
- ・激凹みから立ち直る方法
- ・1つだけ過去を変えられるとしたら?
- ・【あるあるbot連動企画】あるあるbotに投稿したけど採用されなかったあるある募集
- ・【あるあるbot連動企画】フォロワー20万人のアカウントであなたのあるあるを披露してみませんか?
- ・映画のエンドロール観る派?観ない派?
- ・海外旅行から帰ってきたら、まず何を食べる?
- ・誕生日にもらった意外なもの
- ・天使と悪魔選手権
- ・ちょっと先の未来クイズ第2問
- ・【大喜利】【投稿~9/7】 ロボットの住む世界で流行ってる罰ゲームとは?
- ・推しミネラルウォーターはありますか?
- ・都道府県穴埋めゲーム
- ・この人頭いいなと思ったエピソード
- ・準・究極の選択
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
グラフの交点の求め方(Excel)
-
マインクラフト(pc版)で座標...
-
3次元空間上の2つの座標から...
-
エクセルである点からの距離で...
-
2D座標を3D座標に変換する...
-
エクセルで回転する座標の出し方
-
C++でコマンドプロンプトに図形...
-
閉図形の座標の配列が右回りか...
-
最小二乗平面
-
空間上の二点を結ぶ直線上に任...
-
交差する2線分の交点座標の求め方
-
ExcelやAccessで社内の端末の配...
-
凸型の多角形の座標
-
N88-BASICのグラフィック、図形...
-
プログラム
-
WM_NCHITTESTの流れ
-
エクセルシート上のマウスポイ...
-
オートシェイプ円弧の中心点、...
-
多角形の内部かどうか判定する方法
-
Delphiで後ろにあるTPanelや重...
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
グラフの交点の求め方(Excel)
-
マインクラフト(pc版)で座標...
-
3次元空間上の2つの座標から...
-
ダイアログ内コントロールの位...
-
エクセルである点からの距離で...
-
画像ファイルに座標が記録され...
-
始点、終点の二つの座標と半径...
-
直線上にある点の座標の求め方
-
エクセルで回転する座標の出し方
-
以下のプログラムは重心を求め...
-
閉図形の座標の配列が右回りか...
-
シーケンサー(PLC?)で制...
-
座標を持った平面範囲に座標を...
-
円弧の描画について
-
多角形の内部かどうか判定する方法
-
図形が重なりあっているかどうか
-
ガウシアンフィルタのCプログラム
-
ワード上Shapeの位置情報を統一...
-
ピクチャボックスの座標取得
-
交差する2線分の交点座標の求め方
おすすめ情報