三角形の内接円と傍接円の描写

# はじめに
三角形の内心は、三角形の内部にあり、3辺から等距離にある。一方で、傍心は三角形の外部にあり、3辺の延長直線上から等距離にある。したがって、今回は点と直線の距離の公式を用いて、三角形に内接する円と傍接する円を描写する。具体的には、3辺から等距離に点が存在するときに目的関数が最小値0となるような関数を設定して、その目的関数を再急降下法で最適化し、内心や傍心の座標を推定する。ただし、内心は傍心とはことなり、3点存在するので、どれか1つを描写することを目指す。

内接円
![3点を通る内接円の推定.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/c40be4fc-e711-427c-af27-d9669dfbf7a0.png)


傍接円
まず、1つの傍接円を描くと以下のようになる。
![3点を通る傍接円の推定.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/773cc293-c27b-4882-a6fb-96bfca5ad998.png)



次に、下記のように、傍接円は合計で3つ描くことができる。

![3点を通る傍接円の推定.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/7cd6ba97-9524-4c7a-802f-7f526268da41.png)

ちなみに、外接円を描写するプログラムは以下の記事に記載している。

https://qiita.com/arairuca/items/fbecad249f8d4dced499

# 再急降下法による近似解（内心、傍心の座標）の推定
## 導入
3点$A(x_1,y_1),B(x_2,y_2),C(x_3,y_3)$を用いた$\Delta ABC$を考える。

内心もしくは傍心の推定座標を$D(x_0,y_0)$とする。

ところで、直線$AB$の方程式は以下のようにあらわせる。(簡略化のため細かい例外は除外する)

```math
y-y_1=\frac{y_2-y_1}{x_2-x_1}(x-x_1)
```

ゆえに、直線の方程式の一般形は以下のようにあらわせる。

```math
(y_2-y_1)x-(x_2-x_1)y+(x_2y_1-x_1y_2)=0
```

したがって直線$AB$と点$D$の距離$r_1$ は以下のようにあらわせる。

```math
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を取る。

```math
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)$を算出する。


``` math
a_{k+1}=a_{k}-\alpha(\frac{\delta }{\delta a}f(a_{k},b_{k}))
```

``` math
b_{k+1}=b_{k}-\beta(\frac{\delta }{\delta b}f(a_k,b_k))
```
ただし$\alpha,\beta$は$(a,b)$の学習率である。
これは、ボールが坂に下り落ちる様に似ている。偏微分係数は各方向の坂の傾きと捉えることができる。
ただし、今回のプログラムでは簡単のため、$x_0,y_0$を推定するとき、共通の学習率$\alpha$を用いた。

# プログラム
上記の考察をもとに、以下のようなプログラムを作成した。
## 内心
初期位置を三角形の内部に設定することで、三角形の内心の座標の近似値を探索することができる。（一方で、外部に初期位置を設定すれば傍心の座標の近似値を探索することができる。)

```python naisin.py
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()
```

上記のプログラムを実行すると以下のようなグラフが出力される。

### 目的関数の値の推移

![推定内心における目的関数の推移.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/1063a361-89c1-4b9c-83a4-1c0c0bb6163a.png)

このように、試行回数を増やすほど、目的関数は低下していく。
最適化された目的関数による推定座標を用いて、推定内接円を描写すると以下のようになる。

### 内接円の図示

![3点を通る内接円の推定.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/9edf14e1-d492-424b-bb57-b69aac54ed75.png)

## 傍心
次に、初期位置が三角形の外部にある場合について考える。上記のプログラムの初期座標を1セットだけ考えて、それを三角形の外部に配置する。こうすることで、3つある傍心のうちでその初期位置に近い傍心の近似値を算出することができるようになる。

```python bousin.py
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()
```
### 目的関数の値の推移
試行回数は内心のときと比較して増加するが、以下のように目的関数は最適化される。

![推定傍心における目的関数の推移.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/1eb9c6d3-f954-4d3e-96ba-a6e8a1256b37.png)

これは、内心と比較して傍心のほうが探索する範囲は大きくなるからであると考えることができる。
### 傍接円の図示
最適化された目的関数による推定座標を用いて推定傍接円を描写すると以下のようになる。

![3点を通る傍接円の推定.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/9cdfde79-062f-45db-8068-1ae948864b69.png)

### 3つの傍接円を描くプログラム
3つの初期座標を三角形の外部に設定する。ただし、それぞれの3点は同じ辺を最も近い辺としないことにする。

```python bousin3.py
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()
```
### すべての傍接円の図示
これを実行すると以下のようなグラフが出力される。

![3点を通る傍接円の推定.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/7cd6ba97-9524-4c7a-802f-7f526268da41.png)

このように、1つの三角形には3つの傍接円が存在する。

# まとめ

今回は、内心や傍心を推定して、内接円や傍接円を描写するプログラムを作成した。ただし、今回のプログラムは、代数的に求めた厳密解を用いるのではなく、再急降下法による推定値を用いているため、あくまで近似的なものである。というのは、厳密解を用いた場合、外心同様に式が複雑になるからである。ただ、それでも視覚的には十分に精度の良い結果を得ることができた。このように、数値解析や最適化計算といった計算機科学は解析学というツールを合わせて用いると、幾何学とも相性がいい。