四面体の外接球を求める
はじめに
三次元ドロネーを計算する過程において四面体の外接球を計算する必要がありました。（三次元ドロネーに関してはこちら)
結論から言うとWolframMathWorldに記述があるこちらの計算式をPythonで書き下すことで実現しました。
結構しらべてたんですが、ことごとく僕には難解なものばかりでして。。。
WolflramMathWorldに数式がありとても助かったと共に、コード化することに少しは意義ありそうだなと思いましたので記事にまとめます。
理論的なところは抜きで、純粋にコードだけなのでそこは勘弁してください。
実装！
実装はPythonで行っております。
Tetrahedronクラスを作成し、そのクラスメソッドとして外接球の計算を行っています。
Tetrahedronクラスではクラス変数としてPointクラスで作成されたインスタンスを保持する構造を採っています。
実行ファイルでPointクラスとTetrahedronクラスをimportし各インスタンスを作成することで外接球を計算することができます。
今回は実行環境にRhinocerosという3DCADを使用しています。
そちらは今回のメインでは無いのでMain.pyに関する記述は省略させていただきます。。
結果は以下のイメージ。赤い線が四面体です。
Pointクラス
こちらは各座標をクラス変数として保持しているというシンプルな構造です。
# coding: utf-8
class Point:
def __init__(self, point):
self.point = point
self.x = None
self.y = None
self.z = None
self.coordinate = None
def change_coordinate(self):
'''外接円を構成する点を変換するメソッド'''
p = self.point
self.x = p[0]
self.y = p[1]
self.z = p[2]
self.coordinate = [p[0], p[1], p[2]]
Tetrahedronクラス
完全にWolframMathWorldにある数式をPythonコードに落とし込んだだけのものになります。
一点だけ注意事項が有るのですが、僕の使用しています実行環境の問題でnumpyが使用できなかったため、行列計算を行う計算式を自作しています。実際はnumpyを使用したほうが良いです。
# coding: utf-8
import math
class Tetrahedron:
def __init__(self, p1, p2, p3, p4):
self.p1 = p1
self.p2 = p2
self.p3 = p3
self.p4 = p4
self.radius = None
self.center_p = None
def cul_center_p_and_radius(self):
'''外接球の中心点と半径を計算'''
a = [
[self.p1.x, self.p1.y, self.p1.z, 1],
[self.p2.x, self.p2.y, self.p2.z, 1],
[self.p3.x, self.p3.y, self.p3.z, 1],
[self.p4.x, self.p4.y, self.p4.z, 1]
]
d_x = [
[pow(self.p1.x, 2) + pow(self.p1.y, 2) + pow(self.p1.z, 2), self.p1.y, self.p1.z, 1],
[pow(self.p2.x, 2) + pow(self.p2.y, 2) + pow(self.p2.z, 2), self.p2.y, self.p2.z, 1],
[pow(self.p3.x, 2) + pow(self.p3.y, 2) + pow(self.p3.z, 2), self.p3.y, self.p3.z, 1],
[pow(self.p4.x, 2) + pow(self.p4.y, 2) + pow(self.p4.z, 2), self.p4.y, self.p4.z, 1]
]
d_y = [
[pow(self.p1.x, 2) + pow(self.p1.y, 2) + pow(self.p1.z, 2), self.p1.x, self.p1.z, 1],
[pow(self.p2.x, 2) + pow(self.p2.y, 2) + pow(self.p2.z, 2), self.p2.x, self.p2.z, 1],
[pow(self.p3.x, 2) + pow(self.p3.y, 2) + pow(self.p3.z, 2), self.p3.x, self.p3.z, 1],
[pow(self.p4.x, 2) + pow(self.p4.y, 2) + pow(self.p4.z, 2), self.p4.x, self.p4.z, 1]
]
d_z = [
[pow(self.p1.x, 2) + pow(self.p1.y, 2) + pow(self.p1.z, 2), self.p1.x, self.p1.y, 1],
[pow(self.p2.x, 2) + pow(self.p2.y, 2) + pow(self.p2.z, 2), self.p2.x, self.p2.y, 1],
[pow(self.p3.x, 2) + pow(self.p3.y, 2) + pow(self.p3.z, 2), self.p3.x, self.p3.y, 1],
[pow(self.p4.x, 2) + pow(self.p4.y, 2) + pow(self.p4.z, 2), self.p4.x, self.p4.y, 1]
]
a = dit_4(a)
d_x = dit_4(d_x)
d_y = -(dit_4(d_y))
d_z = dit_4(d_z)
center_p_x = d_x / (2 * a)
center_p_y = d_y / (2 * a)
center_p_z = d_z / (2 * a)
self.center_p = [float(center_p_x), float(center_p_y), float(center_p_z)]
distance = math.sqrt(
pow((self.p1.x - self.center_p[0]), 2) + pow((self.p1.y - self.center_p[1]), 2) + pow((self.p1.z - self.center_p[2]), 2))
self.radius = float(distance)
def dit_2(dit):
'''二次元行列の計算'''
cul = (dit[0][0] * dit[1][1]) - (dit[0][1] * dit[1][0])
return cul
def dit_3(dit):
'''三次行列の計算'''
a = [
[dit[1][1], dit[1][2]],
[dit[2][1], dit[2][2]]
]
b = [
[dit[0][1], dit[0][2]],
[dit[2][1], dit[2][2]]
]
c = [
[dit[0][1], dit[0][2]],
[dit[1][1], dit[1][2]]
]
cul = (dit[0][0] * dit_2(a)) - (dit[1][0] * dit_2(b)) + (dit[2][0] * dit_2(c))
return cul
def dit_4(dit):
'''4次元の行列'''
a = [
[dit[1][1], dit[1][2], dit[1][3]],
[dit[2][1], dit[2][2], dit[2][3]],
[dit[3][1], dit[3][2], dit[3][3]]
]
b = [
[dit[0][1], dit[0][2], dit[0][3]],
[dit[2][1], dit[2][2], dit[2][3]],
[dit[3][1], dit[3][2], dit[3][3]]
]
c = [
[dit[0][1], dit[0][2], dit[0][3]],
[dit[1][1], dit[1][2], dit[1][3]],
[dit[3][1], dit[3][2], dit[3][3]]
]
d = [
[dit[0][1], dit[0][2], dit[0][3]],
[dit[1][1], dit[1][2], dit[1][3]],
[dit[2][1], dit[2][2], dit[2][3]]
]
cul = (dit[0][0] * dit_3(a)) - (dit[1][0] * dit_3(b)) + (dit[2][0] * dit_3(c)) - (dit[3][0] * dit_3(d))
return cul
おわりに
最後まで見ていただきありがとうございます。
因みになんですが、四面体をそもそも構成できない位置、つまり同一の直線上に3点が載っている、点が重なっているetc...の条件だとZero Division Error的なモノが返ってきます。ご注意ください。（例外処理書けよって感じですよね笑）