三角形の3つの座標が与えられたとき、ここから内心、外心、垂心、重心および傍心の5心の座標を計算します。
前提
まずは前提として、簡単な関数を定義します。
#1 2つの長さを求める
def length(p1,p2):
s=0
for i,j in zip(p1,p2):
s+=(i-j)**2
return s**0.5
#2 三角形の3つの座標から3辺の長さを求める
def tri(pa,pb,pc):
a=length(pc,pb)
b=length(pa,pc)
c=length(pa,pb)
return a,b,c
#3 三角形の3つの座標からヘロンの公式を使い面積を求める
def square(a,b,c):
s=(a+b+c)/2
return (s*(s-a)*(s-b)*(s-c))**0.5
#4 三角形の3つの長さから、1番目の引数の角のsinの値を求める
def sin_l(a,b,c):
s=square(a,b,c)
return(s*2/b/c)
#5 三角形の3つの長さから、1番目の引数の角のcosの値を求める
def cos_l(a,b,c):
return(b**2+c**2-a**2)/(2*b*c)
#6 原点から2つの点となす角のsinの値を求める
def sin_v(p1,p2):
a,b,c=tri((0,0),p1,p2)
return sin_l(a,b,c)
#7 原点から2つの点となす角のcosの値を求める
def cos_v(p1,p2):
a,b,c=tri((0,0),p1,p2)
return cos_l(a,b,c)
#8 三角形の3つの座標から3つの角のsinの値を求める
def sin_t(pa,pb,pc):
a,b,c=tri(pa,pb,pc)
s=square(a,b,c)
return s*2/b/c,s*2/c/a,s*2/a/b
#9 三角形の3つの座標から3つの角のcosの値を求める
def cos_t(pa,pb,pc):
a,b,c=tri(pa,pb,pc)
return cos_l(a,b,c),cos_l(b,c,a),cos_l(c,a,b)
次に5心を求める上で、強力な関数を定義します。この関数は三角形の3つの点A(xa,ya),B(xb,yb),C(xc,yc)と、三角形の中の点Pに対して△BCP、△CAP、△ABPの面積がga:gb:gcの比率であるとき、点Pの座標は次の通りになります。
$P=\displaystyle \left(\frac{g_{a} x_{a}+ g_{b}x_{b}+ g_{c}x_{c}}{g_{a}+g_{b}+ g_{c}}, \ \frac{g_{a} y_{a}+ g_{b}y_{b}+ g_{c}y_{c}}{g_{a}+g_{b}+ g_{c}}\right)$
def center(pa,pb,pc,ga,gb,gc):
ct=[]
for n in pa:
ct.append(n*ga)
for i,n in enumerate(pb):
ct[i]+=n*gb
for i,n in enumerate(pc):
ct[i]+=n*gc
c=[]
for k in ct:
c.append(k/(ga+gb+gc))
return c少しあか抜けないプログラムで何とかならないかとも思いますが、この関数を使うと内心や外心などを求めるときに、それらが三角形の面積をどの比率で分割されるかを考えればよいわけです。
5心の計算
内心 (incenter)
内心は内接円の中心になり、各辺に下した垂線の長さが内接円の半径になります。内接円の半径は3つに分割した三角形の高さで等しくなるので、面積の比は各辺の長さになります。このため、座標からtri関数を使い、3辺の長さを計算したうえでcenter関数を使います。
def incenter(pa,pb,pc):
a,b,c=tri(pa,pb,pc)
return center(pa,pb,pc,a,b,c)外心 (circumcenter)
外心は、三角形の面積をsin(2A):sin(2B):sin(2C)に分割します。ところで、$sin(2\theta)=2sin(\theta)cos(\theta)$なので、sin_l関数とcos_l関数を利用します。
def circumcenter(pa,pb,pc):
a,b,c=tri(pa,pb,pc)
return center(pa,pb,pc,
sin_l(a,b,c)*cos_l(a,b,c),
sin_l(b,c,a)*cos_l(b,c,a),
sin_l(c,a,b)*cos_l(c,a,b))垂心 (orthocenter)
垂心は、三角形の面積をtan(A):tan(B):tan(C)に分割します。ところで、$tan(\theta)=\frac{sin(\theta)}{cos(\theta)}$なので、sin_l関数とcos_l関数を利用します。
def orthocenter(pa,pb,pc):
a,b,c=tri(pa,pb,pc)
return center(pa,pb,pc,
sin_l(a,b,c)/cos_l(a,b,c),
sin_l(b,c,a)/cos_l(b,c,a),
sin_l(c,a,b)/cos_l(c,a,b))tanに関する関数を定義するのが面倒くさいので、sinとcosを使ったら、外心との比較でcosを掛けるか割るかの違いであるという非常に興味深いことが分かります。
重心 (centergravity)
重心は文字通り、3角形を同じ面積比(1:1:1)で分割するので、とても単純です。
def centergravity(pa,pb,pc):
a,b,c=tri(pa,pb,pc)
return center(pa,pb,pc,1,1,1)傍心
傍心は、面倒くさそうですが、内心の考え方を延長すると計算することができます。傍心は3つあるのでリストを3つまとめて返すようにします。
def excenter(pa,pb,pc):
a,b,c=tri(pa,pb,pc)
return (center(pa,pb,pc,-a,b,c),
center(pa,pb,pc,a,-b,c),
center(pa,pb,pc,a,b,-c))実行
点A(12,18),B(390,18),C(102,138)からなる三角形で、実行してみました。
x=(12,18),(390,18),(102,138)
print('内心',incenter(*x))
print('外心',circumcenter(*x))
print('垂心',orthocenter(*x))
print('重心',centergravity(*x))
print('傍心',excenter(*x)) 内心 [120.0, 72.0]
外心 [201.00000000000003, -30.0]
垂心 [102.0, 234.0]
重心 [168.0, 58.0]
傍心 ([432.0, 228.0], [-30.0, 102.0], [282.0, -522.0])実行結果については、sympyモジュールで一部確認することができます。
import sympy.geometry
A, B, C = sympy.geometry.Point(12, 18), sympy.geometry.Point(390, 18), sympy.geometry.Point(102, 138)
tri = sympy.geometry.Triangle(A, B, C)
print('内心 = ',tri.incenter)
print('半径 = ',tri.inradius)
print('外心 = ',tri.circumcenter)
print('垂心 = ',tri.orthocenter)内心 = Point2D(120, 72)
半径 = 54
外心 = Point2D(201, -30)
垂心 = Point2D(102, 234)内心、外心、垂心については、結果が正しいことが確認できました。重心は簡単すぎるのか計算方法が見つからず、傍心についてはexcenterというプロパティがありますが、目的が違うようです。ただ、Sympyには内心の半径を求めるなど異なった機能もあるので、この辺は別に考えます。
