2直線の最短距離となる座標を算出したいです。

2直線の最短距離座標を図を参考に算出しました。

係数s、tはsymsを用いて定義しました。

しかし、データ数が627912個もあります。

matlabを走らせると計算時間が12時間ほどかかりました。

もっと簡単な算出方法、短時間で計算ができるようなスクリプトを作りたいのですが、中々上手くいきません。

良い方法があれば教えてほしいです。

Answer 1

[Kazuya]

627912点のデータか何のデータか分かりませんが、2直線の式をまず求めて

ねじれの位置にある2直線間の最短距離

で距離を計算するのが計算負荷が小さそうですがどうでしょ。

[ryo tanaka]

回答ありがとうございます。

627912個とは、最短距離を求める回数でもあります。

左右の2直線は、左が972個の直線があり、右は646個の直線があり、左右の各ペアどうしの差に短距離を算出したいので、972×646＝627912個です。

資料を添付していただきありがとうございます。

私もこの添付資料を参考に計算式を作成しましたが、係数を用いているので、計算時間が長くなってしまいます。

[Kazuya]

なるほど、627912 個の組み合わせについて計算するんですね。

もし 5 組程度の組み合わせデータを見せてもらえれば、処理をベクトル化するなどできるかどうか考えてみますけどいかがでしょう？

[ryo tanaka]

ありがとうございます。データを添付します。

2直線の組合せはAO直線とBC直線です。O点とC点は固定の座標です。

% B,O,C,A = [x座標,y座標,z座標] 
B=[39.8125,-0.9625,-198;39.9875,-1.8375,-198;39.9875,-0.9625,-198;39.9875,-0.7875,-198;40.1625,-7.2625,-198];
O=[-41.2125,-0.4375,-198;-41.2125,-0.2625,-198;-41.0375,-6.2125,-198;-41.0375,-0.9625,-198;-41.0375,-0.7875,-198];
C=[25,0,0];
A=[-25,0,0];

[Yoshio]

AとＯが入れ替わっていて、O = [25, 0, 0]が正しいですか？

[Kazuya]

画像も見えて状況がつかめてきました。Yoshio さんのいう通り matlabFunction を使う方法がよいのかも？

係数s、t を syms で定義して処理したコードも見せて頂くことは可能ですか？考えるベースにできればと思いまして。例えば上で頂いたデータに対して適用した形とか。。

[ryo tanaka]

回答ありがとうございます。

返信が遅くなり、大変申し訳ありませんでした。

係数s、t を syms で定義して処理したコードは添付画像通りに算出していますので、、

コードを添付した方がよろしいでしょうか？

Answer 2

[Yoshio]

もし、最短距離だけを求めるのであれば、

2直線の距離 - 高精度計算サイト

から、,, とすればベクトル化した計算式ができそうです。

[Yoshio]

とりあえずC=[25,0,0];として距離だけ求めるコードを書いてみました。

ベクトル化のために、データは縦に[x;y;z]と並べ直しています。

ただし、この例はBとＡの各列ごとに直線が対応しているとしてのみ計算しています(Bi列にはＡi列のみ計算）。

B=[39.8125,-0.9625,-198;39.9875,-1.8375,-198;39.9875,-0.9625,-198;39.9875,-0.7875,-198;40.1625,-7.2625,-198];
A=[-41.2125,-0.4375,-198;-41.2125,-0.2625,-198;-41.0375,-6.2125,-198;-41.0375,-0.9625,-198;-41.0375,-0.7875,-198];
C=[25,0,0];
O=[-25,0,0];

%%
A = A';
B = B';
C = C';
O = O';

tic
v1 = A-O;
v2 = B-C;
P12 = O-C;
v1xv2 = cross(v1,v2);
d = abs(P12'*v1xv2./vecnorm(v1xv2));
toc

後は、BとＡの組み合わせのを作成して同様にすれば良さそうですが、とりあえず

n = 1000000;
A = rand(3,n);
B = rand(3,n);

として計算したところ2.3 GHz Intel Core i5のMac Miniでも0.13秒で終わりました

[Yoshio]

AとBに各々独立な直線を指定する点を入れるようにして、直線の組み合わせも考慮すると

C=[25,0,0];
O=[-25,0,0];
tic
A = rand(3,972);
B = rand(3,646);

AA = repmat(A,1,size(B,2));
BB = reshape(repmat(B,size(A,2),1),3,size(B,2)*size(A,2));

C = C';
O = O';

v1 = AA-O;
v2 = BB-C;
P12 = O-C;
v1xv2 = cross(v1,v2);
d = abs(P12'*v1xv2./vecnorm(v1xv2));
toc

972×646＝627912の組だと0.1秒でした。

[Kazuya]

ありがとうございます。勉強になりました！

[ryo tanaka]

回答していただきありがとうございます。

返信が遅くなり大変申し訳ありませんでした。

AとOが入れ替わっていました。すみません。

O = [25, 0, 0]が正しいです。

処理時間0.1秒というのは自分の理想をはるかに超えるはやさですね。

一度このコードを用いて計算させてみたいと思います。

最終的には最短距離と2直線の最短距離となる座標も求めたいと思っています。。

[Yoshio]

確認ありがとうございます。O = [25, 0, 0]で良かったです。

何回か投稿を編集されているようで、私が見たときは最短距離だけ出ればよいかと思ったのですが、座標も必要なのですね。

その場合、t,s を求める計算をv1, v2, P12等で記述できれば高速化可能と思います。

[ryo tanaka]

何度か投稿の編集を行っていましたので、伝えきれていなかったと思います。

失礼いたしました。

なるほど、t、sの記述方法を変えることで計算時間の短縮ができそうですね。

Answer 3

[Yoshio]

添付の図が見えませんので、良く分からない所がありますが、シンボリックが解をmatlabFunctionに変換できるので、一度数式で得られた解を数値計算向けの関数にして利用するのはどうでしょうか？　一度試してみてください。

[ryo tanaka]

回答ありがとうございます。添付した図が見れず申し訳ありません。修正しました。

matlabFunctionは知りませんでした。一度試してみようと思います。ありがとうございます。

Answer 4

[Yoshio]

直線モデルを以下のように置く

とすると、最短距離の条件から

よって

これをについて解くと

ここで

上記に基づき、ベクトル化を考慮してプログラムすると

% B,O,C,A = [x座標,y座標,z座標] 
B=[39.8125,-0.9625,-198;39.9875,-1.8375,-198;39.9875,-0.9625,-198;39.9875,-0.7875,-198;40.1625,-7.2625,-198];
A=[-41.2125,-0.4375,-198;-41.2125,-0.2625,-198;-41.0375,-6.2125,-198;-41.0375,-0.9625,-198;-41.0375,-0.7875,-198];
C=[25,0,0];
O=[-25,0,0];

A = A';
B = B';
C = C';
O = O';

% A = rand(3,972);
% B = rand(3,646);
AA = repmat(A,1,size(B,2));
BB = reshape(repmat(B,size(A,2),1),3,size(B,2)*size(A,2));

v1 = AA-O;
v2 = BB-C;
P12 = O-C;
v1xv2 = cross(v1,v2);
d = abs(P12'*v1xv2./vecnorm(v1xv2));

tic
x1 = O;
x2 = C;
x12 = x1-x2;
v1v1 = sum(v1.*v1);
v1v2 = sum(v1.*v2);
v2v2 = sum(v2.*v2);
D = -v1v1.*v2v2 + v1v2.^2;
b = -[sum(v1.*x12);sum(v2.*x12)];
t1 = (-v2v2.*b(1,:)+v1v2.*b(2,:))./D;
t2 = (-v1v2.*b(1,:)+v1v1.*b(2,:))./D;
P1 = x1+t1.*v1;
P2 = x2+t2.*v2;
d0 = vecnorm(P1-P2);
toc
max(abs(d0-d))
min(abs(D))

972×646＝627912の組だと約0.2秒でした

dとd0の差はepsの10倍程なので、合っているかと思います。なお、2つの直線が全く平行の場合は、D=0となり、解は不定となるので、この場合は別途考慮する必要がありますね。

[ryo tanaka]

とても分かりやすい回答ありがとうございます。

このコードだと2点座標及び最短距離が求められるという事ですね。

自分のプログラムに用いて計算させてみたいと思います。

勉強不足で申し訳ありません、一つわからないことがあるのですが、、

d0を求める際に用いているvecnormがよくわかりません。

vecnormの解説を拝見しましたが、あまり理解ができませんでした。

もしよろしければ教えて頂けないでしょうか？

[Yoshio]

今回の場合、vecnormの解説にありますが

”N = vecnorm(A) は、A の 2 ノルム、つまりユークリッド ノルムを返します。

• A が行列の場合、vecnorm は各列のノルムを返します。

を利用しています。（P1-P2）は2つの座標の差ベクトル（3次元）を627912個並べた3x627912の行列になっているので、vecnorm(P1-P2)を計算することで、各列ベクトルのノルム（＝2点P1とP2の距離）が627912個並ぶことになります。

まず適当な行列を作って、vecnormを試してみることをお勧めします。

[ryo tanaka]

返信遅くなり申し訳ありません。

各列ベクトルの2点間距離を求めることができるという事ですね。

分かりやすい説明ありがとうございました。

一度試してみたいと思います。

[ryo tanaka]

自分のデータを紹介して頂いたコードで走らせてみました。

計算時間も早く、P1,P2,ｄ(2点間距離)の算出ができました。

ありがとうございました。

また、一つどうしても分からないコードがあったのですが、

下記のコードの4,5行目の計算式がどういうことなのかわかりません。。

ｄの結果自体はあってそうなので、式は正しいのだろうと理解しているのですが、、

もしよろしければ教えて頂きたいです。

v1 = AA-O;
v2 = BB-C;
P12 = O-C;
v1xv2 = cross(v1,v2);
d = abs(P12'*v1xv2./vecnorm(v1xv2));

[Yoshio]

結果のご報告ありがとうございます。うまく動いたようでよかったです。

コードの4、5行目の式は、最初回答で示しました2直線の距離 - 高精度計算サイトの以下の式を参照して作成しました。

「２直線の傾きの外積ベクトルの向きは2つの直線に垂直になります。

この向きの単位ベクトルと各直線上の任意点間のベクトルとの内積の大きさが２直線間の最短距離になります。」

[ryo tanaka]

ご丁寧に分かりやすい説明をして頂いて感謝いたします。

ようやく理解できました。ありがとうございます。

計算で求めたｄは、

[1点目のAO直線に対するBC直線1～end、2点目のAO直線に対するBC直線1～end、・・・、end点目のAO直線に対するBC直線1～end]

という順番でならんでいるということであっていますでしょうか？

質問が多くて申し訳ありません、、、

[Yoshio]

[1点目のAO直線に対するBC直線1～end、2点目のAO直線に対するBC直線1～end、・・・、end点目のAO直線に対するBC直線1～end]

という順番でならんでいる

という記述の意味がよくわからないのですが、MATLABの良さは、アイディアを簡単に試せる点にありますので、サンプルデータを作って、確認していただくのが良いかと思います。

[ryo tanaka]

質問の仕方が悪かったです。すみません。

簡単なサンプルデータで確認してみます。