【Python】ねじれの位置にある２直線の最短距離を計算する方法

こんにちは、ヒガシです。

このページではPythonを使って、ねじれの位置にある２本の直線の最短距離（共通垂線の長さ）を計算する方法をご紹介していきます。

基本的にあなたがやることは上の画像に示すようにそれぞれの線を構成する２点の座標（a,b,c,d）を指定するだけでOKです。

それではさっそくやっていきましょう！

ねじれの位置にある２直線の最短距離の計算方法
Pythonで共通垂線の長さを計算するサンプルコード
サンプルコードの実行結果確認
プログラム実行上の注意点
おわりに

ねじれの位置にある２直線の最短距離の計算方法

ねじれの位置にある２直線の最短距離を計算する方法については以下のサイトが非常に参考になりましたので掲載しておきます。

空間の2直線の最短距離

上野竜生です。ねじれの位置にある空間上の2直線上にそれぞれ点P,QをとったときのPQの最小を考えましょう。裏…

今回は計算結果が正解かどうか確認するために上記参考サイトと全く同じ２本の線に対して計算を実施していきます。

なお、数学的に計算する場合は上記サイトのように計算していけばいいですが、コンピューターがこの計算をするのはあまり得意ではありません。

というわけここからは、コンピュータ上で今回やりたい計算を実施する方法をざっくり解説していきます。

イメージは以下の通りです。

①まずは適当にそれぞれの直線上に点を置き、両者の距離を計測します。

②それぞれの点をすこしだけ動かしてまた距離を計算します。

これを何度も繰り替えし、最短距離を探していく作業をおこなっていきます。

まぁ難しいことは考えず、いったんコーディングしてみましょう。

Pythonで共通垂線の長さを計算するサンプルコード

それではさっそくですが、サンプルコードのご紹介です。

#ライブラリインポート
import numpy as np
import matplotlib.pyplot as plt 
from mpl_toolkits.mplot3d import Axes3D

#ab,cdでそれぞれ繋がる２直線の最短距離を計算する
#参照(https://math-juken.com/kijutu/kyoutusuisen/)
a=[0,0,0]
b=[1,2,3]
c=[4,5,-6]
d=[7,8,9]
line1=[a,b]
line2=[c,d]

#２直線の距離を計算する関数
def calc_distance(s,t):
    a,b=line1[0],line1[1]
    c,d=line2[0],line2[1]
    ABbec=np.array([b[0]-a[0],b[1]-a[1],b[2]-a[2]])
    CDbec=np.array([d[0]-c[0],d[1]-c[1],d[2]-c[2]])
    PQbec=(c+CDbec*s)-(a+ABbec*t)
    distance=np.linalg.norm(PQbec)
    return distance

#２直線の距離を探索する関数
def make_map(ss,tt):
    z_list=[]
    ij_list=[]
    for i in range(len(ss)):
        for j in range(len(tt)):
            s,t=ss[i],tt[j]
            z=calc_distance(s,t)
            z_list.append(z)
            ij_list.append([i,j])
    return z_list,ij_list

#大きな範囲で最小値をざっくり探索
ini_max=2000
ini_min=-2000
delta=200
ss=np.linspace(ini_min,ini_max,delta)
tt=np.linspace(ini_min,ini_max,delta)
z,ij=make_map(ss,tt)
min_val=np.amin(np.array(z))
min_index=z.index(min_val)

#範囲を狭めながら最小値を繰り返し計算
minz_list=[]
for k in range(7):
    i,j=ij[min_index][0],ij[min_index][1]
    smax,smin=min(ss[i+2],ini_max),max(ss[i-2],ini_min)
    tmax,tmin=min(tt[j+2],ini_max),max(tt[j-2],ini_min)
    ss=np.linspace(smin,smax,delta)
    tt=np.linspace(tmin,tmax,delta)
    z,ij=make_map(ss,tt)
    min_val=np.amin(np.array(z))
    min_index=z.index(min_val)
    minz_list.append(min_val)
    print('min_distance=' + str(min_val))

#正解データの計算(https://math-juken.com/kijutu/kyoutusuisen/)
P=np.array([5,10,15])
Q=np.array([73,82,131])/9
PQ=P-Q
print('answer=' +str(np.linalg.norm(PQ)))

非常に長いですが、基本的にあなたがやることは８～１１行目のa,b,c,d各点の座標を入力するだけでOKです。