はじめに
三角形の内心は、三角形の内部にあり、3辺から等距離にある。一方で、傍心は三角形の外部にあり、3辺の延長直線上から等距離にある。したがって、今回は点と直線の距離の公式を用いて、三角形に内接する円と傍接する円を描写する。具体的には、3辺から等距離に点が存在するときに目的関数が最小値0となるような関数を設定して、その目的関数を再急降下法で最適化し、内心や傍心の座標を推定する。ただし、内心は傍心とはことなり、3点存在するので、どれか1つを描写することを目指す。
次に、下記のように、傍接円は合計で3つ描くことができる。
ちなみに、外接円を描写するプログラムは以下の記事に記載している。
再急降下法による近似解(内心、傍心の座標)の推定
導入
3点$A(x_1,y_1),B(x_2,y_2),C(x_3,y_3)$を用いた$\Delta ABC$を考える。
内心もしくは傍心の推定座標を$D(x_0,y_0)$とする。
ところで、直線$AB$の方程式は以下のようにあらわせる。(簡略化のため細かい例外は除外する)
y-y_1=\frac{y_2-y_1}{x_2-x_1}(x-x_1)
ゆえに、直線の方程式の一般形は以下のようにあらわせる。
(y_2-y_1)x-(x_2-x_1)y+(x_2y_1-x_1y_2)=0
したがって直線$AB$と点$D$の距離$r_1$ は以下のようにあらわせる。
r_1=\frac{|(y_2-y_1)x_0-(x_2-x_1)y_0+(x_2y_1-x_1y_2)|}{\sqrt{(y_2-y_1)^2+(x_2-x_1)}}
同様に、直線$BC$と点$D$の距離を$r_2$,直線$AC$と点$D$の距離を$r_3$とする。もし、点$D$が内心もしくは傍心の場合は、以下に定義する目的関数が最小値0を取る。
f(x_0,y_0)=(r_1-r_2)^2+(r_2-r_3)^2+(r_3-r_1)^2
したがって上記の目的関数$f(x_0,y_0)$を最適化するような座標を初期値を与えて、最急降下法で推定する。(初期値が三角形の内部にあれば内心になりやすく、外部にあれば傍心になりやすくなる。)
最急降下法
ここでも簡潔に最急降下法アルゴリズムを説明しておこう。
$z=f(x,y)$の最小値を求めることを考える。
解空間の最小値を取る極値が一つ解空間の谷が一つの場合に有効な最小値の推定手法である。$(a,b)$に初期値を与え、以下のような式を用いて$(a,b)$を更新していく。つまり、$(a_k,b_k)$を算出する。
a_{k+1}=a_{k}-\alpha(\frac{\delta }{\delta a_{i}}f(a_{k},b_{k}))
b_{k+1}=b_{k}-\beta(\frac{\delta }{\delta b_{k}}f(a_k,b_k))
ただし$\alpha,\beta$は$(a,b)$の学習率である。
これは、ボールが坂に下り落ちる様に似ている。偏微分係数は各方向の坂の傾きと捉えることができる。
ただし、今回のプログラムでは簡単のため、$x_0,y_0$を推定するとき、共通の学習率$\alpha$を用いた。
プログラム
上記の考察をもとに、以下のようなプログラムを作成した。
内心
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
#再急降下法の学習率
alpha=0.1
#微小分
delta=1.0*10**-5
#三角形の座標
x1=1
y1=2
x2=-3
y2=4
x3=-5
y3=-6
#内心、傍心の推定値の初期値
x0=0
y0=1
#目的関数の定義
def f(x0,y0):
#内心、傍心の推定値を用いて推定半径を計算
r1=abs((y2-y1)*x0-(x2-x1)*y0+(x2*y1-x1*y2))/(((x2-x1)**2)+(y2-y1)**2)**0.5
r2=abs((y3-y2)*x0-(x3-x2)*y0+(x3*y2-x2*y3))/(((x3-x2)**2)+(y3-y2)**2)**0.5
r3=abs((y1-y3)*x0-(x1-x3)*y0+(x1*y3-x3*y1))/(((x1-x3)**2)+(y1-y3)**2)**0.5
#目的関数の値を計算
X=(r1-r2)**2+(r1-r3)**2+(r2-r3)**2
return X
i=0
f_ary=[]
i_ary=[]
while f(x0,y0)>delta:
#微分の計算
dfdx=(f(x0+delta,y0)-f(x0,y0))/delta
dfdy=(f(x0,y0+delta)-f(x0,y0))/delta
#再急降下法によるパラメータの更新
x0=x0-alpha*dfdx
y0=y0-alpha*dfdy
#反復回数と目的関数の値を格納
f_ary.append(f(x0,y0))
i_ary.append(i)
i=i+1
#結果の出力
plt.xlabel('反復回数')
plt.ylabel('目的関数')
plt.plot(i_ary,f_ary)
#plt.savefig('推定傍心における目的関数の推移.png')
plt.savefig('推定内心における目的関数の推移.png')
plt.show()
#半径の推定
r0=abs((y2-y1)*x0-(x2-x1)*y0+(x2*y1-x1*y2))/(((x2-x1)**2)+(y2-y1)**2)**0.5
#円の描画(極座標)
theta=np.linspace(0,math.pi*2,100)
#パラメータ表示
x=x0+r0*np.cos(theta)
y=y0+r0*np.sin(theta)
#3点の座標をプロット
plt.plot(x1,y1,color='red',marker='o')
plt.plot(x2,y2,color='red',marker='o')
plt.plot(x3,y3,color='red',marker='o')
#plt.plot([x1,x2,x3,x1],[y1,y2,y3,y1],color='red')
plt.axline((x1,y1),(x2,y2),color='red')
plt.axline((x2,y2),(x3,y3),color='red')
plt.axline((x3,y3),(x1,y1),color='red')
#推定した円をプロット
plt.plot(x,y,color='blue')
#アスペクト比を等しく
plt.axis('equal')
#plt.savefig('3点を通る傍接円の推定.png')
plt.savefig('3点を通る内接円の推定.png')
plt.show()
上記のプログラムを実行すると以下のようなグラフが出力される。
このように、試行回数を増やすほど、目的関数は低下していく。
最適化された目的関数による推定座標を用いて、推定内接円を描写すると以下のようになる。
傍心
次に、初期位置が三角形の外部にある場合について考える。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
#再急降下法の学習率
alpha=0.1
#微小分
delta=1.0*10**-5
#三角形の座標
x1=1
y1=2
x2=-3
y2=4
x3=-5
y3=-6
#内心、傍心の推定値の初期値
x0=-5
y0=1
#目的関数の定義
def f(x0,y0):
#内心、傍心の推定値を用いて推定半径を計算
r1=abs((y2-y1)*x0-(x2-x1)*y0+(x2*y1-x1*y2))/(((x2-x1)**2)+(y2-y1)**2)**0.5
r2=abs((y3-y2)*x0-(x3-x2)*y0+(x3*y2-x2*y3))/(((x3-x2)**2)+(y3-y2)**2)**0.5
r3=abs((y1-y3)*x0-(x1-x3)*y0+(x1*y3-x3*y1))/(((x1-x3)**2)+(y1-y3)**2)**0.5
#目的関数の値を計算
X=(r1-r2)**2+(r1-r3)**2+(r2-r3)**2
return X
i=0
f_ary=[]
i_ary=[]
while f(x0,y0)>delta:
#微分の計算
dfdx=(f(x0+delta,y0)-f(x0,y0))/delta
dfdy=(f(x0,y0+delta)-f(x0,y0))/delta
#再急降下法によるパラメータの更新
x0=x0-alpha*dfdx
y0=y0-alpha*dfdy
#反復回数と目的関数の値を格納
f_ary.append(f(x0,y0))
i_ary.append(i)
i=i+1
#結果の出力
plt.xlabel('反復回数')
plt.ylabel('目的関数')
plt.plot(i_ary,f_ary)
plt.savefig('推定傍心における目的関数の推移.png')
#plt.savefig('推定内心における目的関数の推移.png')
plt.show()
#半径の推定
r0=abs((y2-y1)*x0-(x2-x1)*y0+(x2*y1-x1*y2))/(((x2-x1)**2)+(y2-y1)**2)**0.5
#円の描画(極座標)
theta=np.linspace(0,math.pi*2,100)
#パラメータ表示
x=x0+r0*np.cos(theta)
y=y0+r0*np.sin(theta)
#3点の座標をプロット
plt.plot(x1,y1,color='red',marker='o')
plt.plot(x2,y2,color='red',marker='o')
plt.plot(x3,y3,color='red',marker='o')
#plt.plot([x1,x2,x3,x1],[y1,y2,y3,y1],color='red')
plt.axline((x1,y1),(x2,y2),color='red')
plt.axline((x2,y2),(x3,y3),color='red')
plt.axline((x3,y3),(x1,y1),color='red')
#推定した円をプロット
plt.plot(x,y,color='blue')
#アスペクト比を等しく
plt.axis('equal')
plt.savefig('3点を通る傍接円の推定.png')
#plt.savefig('3点を通る内接円の推定.png')
plt.show()
試行回数は内心のときと比較して増加するが、以下のように目的関数は最適化される。
最適化された目的関数による推定座標を用いて推定傍接円を描写すると以下のようになる。
3つの傍接円を描くプログラム
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
#再急降下法の学習率
alpha=0.1
#微小分
delta=1.0*10**-5
#三角形の座標
x1=1
y1=2
x2=-3
y2=4
x3=-5
y3=-6
#内心、傍心の推定値の初期値
#x0=-5
#y0=1
#目的関数の定義
def f(x0,y0):
#内心、傍心の推定値を用いて推定半径を計算
r1=abs((y2-y1)*x0-(x2-x1)*y0+(x2*y1-x1*y2))/(((x2-x1)**2)+(y2-y1)**2)**0.5
r2=abs((y3-y2)*x0-(x3-x2)*y0+(x3*y2-x2*y3))/(((x3-x2)**2)+(y3-y2)**2)**0.5
r3=abs((y1-y3)*x0-(x1-x3)*y0+(x1*y3-x3*y1))/(((x1-x3)**2)+(y1-y3)**2)**0.5
#目的関数の値を計算
X=(r1-r2)**2+(r1-r3)**2+(r2-r3)**2
return X
def outer(x0,y0):
i=0
f_ary=[]
i_ary=[]
while f(x0,y0)>delta:
#微分の計算
dfdx=(f(x0+delta,y0)-f(x0,y0))/delta
dfdy=(f(x0,y0+delta)-f(x0,y0))/delta
#再急降下法によるパラメータの更新
x0=x0-alpha*dfdx
y0=y0-alpha*dfdy
#反復回数と目的関数の値を格納
f_ary.append(f(x0,y0))
i_ary.append(i)
i=i+1
#結果の出力
# plt.xlabel('反復回数')
# plt.ylabel('目的関数')
# plt.plot(i_ary,f_ary)
# plt.savefig('推定傍心における目的関数の推移.png')
# #plt.savefig('推定内心における目的関数の推移.png')
# plt.show()
#半径の推定
r0=abs((y2-y1)*x0-(x2-x1)*y0+(x2*y1-x1*y2))/(((x2-x1)**2)+(y2-y1)**2)**0.5
#円の描画(極座標)
theta=np.linspace(0,math.pi*2,100)
#パラメータ表示
x=x0+r0*np.cos(theta)
y=y0+r0*np.sin(theta)
#推定した円をプロット
plt.plot(x,y,color='blue')
#3点の座標をプロット
plt.plot(x1,y1,color='red',marker='o')
plt.plot(x2,y2,color='red',marker='o')
plt.plot(x3,y3,color='red',marker='o')
#plt.plot([x1,x2,x3,x1],[y1,y2,y3,y1],color='red')
plt.axline((x1,y1),(x2,y2),color='red')
plt.axline((x2,y2),(x3,y3),color='red')
plt.axline((x3,y3),(x1,y1),color='red')
outer(-5,-5)
outer(10,10)
outer(0,-5)
#アスペクト比を等しく
plt.axis('equal')
plt.savefig('3点を通る傍接円の推定.png')
#plt.savefig('3点を通る内接円の推定.png')
plt.show()
これを実行すると以下のようなグラフが出力される。
このように、1つの三角形には3つの傍接円が存在する。
まとめ
今回は、内心や傍心を推定して、内接円や傍接円を描写するプログラムを作成した。ただし、今回のプログラムは、代数的に求めた厳密解を用いるのではなく、再急降下法による推定値を用いているため、あくまで近似的なものである。ただ、それでも視覚的には十分に精度の良い結果を得ることができた。