はじめに
前回、任意の座標点が多角形の内側か外側か判定する記事を書きました。
実際は、3軸平面において、多面体の内側か外側かを判定する必要があったので、
2次元から3次元に拡張してみました。
環境
python-3.7.4
numpy-1.17.1
matplotlib-3.1.1
考え方
2次元で作成した方法を少し応用しました。
例
以下の図のように、今回は四面体があるとして、青線で囲っています。（赤が任意の点、四面体の内側に存在する）
赤い点は任意の点として、この点が四面体の内側 or 外側どちらに存在するかを判定します。
実装
まず、モジュールのインポート
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from numpy import linalg as LA
次に、青色で囲った四面体座標を作成します。
後々、平面の重心位置も必要となるので、ついでに出力します。
# 平面内の3点を指定
#面1
p1 = np.array([8,13,19])
p2 = np.array([18,1,4])
p3 = np.array([1,7,1])
# 3点の重心
g1 = (p1 + p2 + p3) / 3
#面2
p4 = np.array([8,13,19])
p6 = np.array([18,1,4])
p5 = np.array([10,18,6])
# 3点の重心
g2 = (p4 + p5 + p6) / 3
#面3
p7 = np.array([18,1,4])
p9 = np.array([1,7,1])
p8 = np.array([10,18,6])
# 3点の重心
g3 = (p7 + p8 + p9) / 3
#面4
p10 = np.array([8,13,19])
p11 = np.array([1,7,1])
p12 = np.array([10,18,6])
# 3点の重心
g4 = (p10 + p11 + p12) / 3
これらの座標点を囲んでいきます。
座標を囲んでいく作業は以下のようなことをしないといけないのかは不明です。
やり方がわからなかったので、下記のような手順にしました。
#座標値の取得
x = []
y = []
z = []
for i in range(1,13):
exec('x_num = p{}[0]'.format(i))
x.append(x_num)
exec('y_num = p{}[1]'.format(i))
y.append(y_num)
exec('z_num = p{}[2]'.format(i))
z.append(z_num)
x.append(p1[0])
y.append(p1[1])
z.append(p1[2])
そして、任意座標を指定して、描画します。
無駄なコードがいくつかある気がします。すみません。
#任意の点を指定
t = np.array([10,15,6])
# FigureとAxes
fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111, projection='3d')
#グラフの描画(一筆描き)
ax.plot(x, y, z, marker="o",markersize=6, color='blue')
#四面体座標の描画
ax.plot(x, y, z, marker="o",markersize=6, color='blue', linestyle = 'none')
#任意座標の描画
ax.scatter(t[0], t[1], t[2],s=40, color='red')
ここまでで、上記の図が得られます。（軸の値は異なります）
ここから、それぞれの面の法線ベクトルを得ます。
以下が簡単な解説です。
ここでは、省略しますが、法線ベクトルの向きは、座標を時計回り or 反時計周りで選ぶ順番によって変わるので、
適宜向かせたい面に向けてください。
①任意座標が内側にある時、法線ベクトルに対して、任意座標と重心を結ぶベクトルの角度は90度以下
②任意座標が外側にある時、90度以上
となります。
角度の求め方は、ベクトルの内積を用います。
ここから、実際のコードになります。
複数平面について、扱うので関数化しています。
#複数平面から算出する角度を入れる箱
list =[]
def point_plane(axes, pa, pb, pc,
xrange, yrange, zrange,
loc = [0, 0, 0],
vcolor="red", pcolor="blue", alpha=0.5):
u = pb - pa
v = pc - pa
#法線ベクトルw
w = np.cross(u, v)
#任意の点へのベクトルh
h = t - loc
#法線ベクトルwと平面から任意点のベクトルhの内積
i = np.inner(w, h)
n = LA.norm(w) * LA.norm(h)
#内積から角度aを求める
c = i / n
a = np.rad2deg(np.arccos(np.clip(c, -1.0, 1.0)))
list.append(a)
次に、作成する平面について関数化しつつ、軸ラベル等定義しておきます。
def nvector_plane(axes, vector, point,
xrange, yrange, zrange,
loc = [0, 0, 0],
vcolor="red", pcolor="blue", alpha=0.5):
# 軸ラベルの設定
axes.set_xlabel("x", fontsize = 16)
axes.set_ylabel("y", fontsize = 16)
axes.set_zlabel("z", fontsize = 16)
# 軸範囲の設定
axes.set_xlim(xrange[0], xrange[1])
axes.set_ylim(yrange[0], yrange[1])
axes.set_zlim(zrange[0], zrange[1])
# 格子点の作成
x = np.arange(xrange[0], xrange[1], 0.2)
y = np.arange(yrange[0], yrange[1], 0.2)
xx, yy = np.meshgrid(x, y)
# 平面の方程式
#point=p1
#vector=w
zz = point[2] - (vector[0]*(xx-point[0])+vector[1]*(yy-point[1])) / vector[2]
必要であれば、下記コードを入れるとプロットする平面や法線ベクトルを可視化できます。
私の場合、見えにくくなるので、削除しました。
ax.plot_surface(xx, yy, zz, color=pcolor, alpha=alpha)
# 法線ベクトルをプロット
axes.quiver(始点x,y,z,w[0],w[1],w[2])
axes.quiver(loc[0], loc[1], loc[2],
vector[0], vector[1], vector[2],
color = vcolor, length = 1, arrow_length_ratio = 0.2)
そして、最後に、作成した関数に座標値を代入していきます。
list内に、角度が入っていきます。
最後に、90度以上 or 以下であるかどうかを判断する処理を加えて終了です。
# 3点を通る平面をプロット
#p1～p3,xr～zr,重心位置gをlocに渡す
point_plane(ax, p1, p2, p3, xr, yr, zr,loc=g1)
point_plane(ax, p4, p5, p6, xr, yr, zr,loc=g2)
point_plane(ax, p7, p8, p9, xr, yr, zr,loc=g3)
point_plane(ax, p10, p11, p12, xr, yr, zr,loc=g4)
print(list)
#リスト内に1つでも90以上があれば、Trueで外側
#なければ、Falseで内側
if any((x >= 90 for x in list)) == True:
print('座標点は立体の外側')
else:
print('座標点は立体の内側')
出力結果は以下の通り。
座標点を変えても対応でき、内外判定が合っていることを確認できました。
総括
今回は、試行錯誤しながら、処理を作っていきました。
四面体で判断できたこと、処理の考え方を明確にしたことで、
多面体においても同様の処理で内外判定ができるだろうと考えます。
今後は実際に、業務に反映して、判断処理の自動化を図っていく予定です。
本記事をご覧いただきありがとうございました。
参考
以下のサイトを参考にさせていただきました。ありがとうございます。
https://python.atelierkobato.com/plane/