sympyとmatplotlibで三角形の五心やオイラー線，九点円などを図示し，フォイエルバッハの定理を確認しよう


この記事は，pythonをかじったことがある人で，matplotlibを初めて使う人（つまり私）を対象に，初等幾何の集大成（は言い過ぎかな？）である九点円の図示および，フォイエルバッハの定理（九点円が内接円と傍接円に接する）が実際に正しいかを目で確認しよう，という趣旨のものです．

SymPyの`sympy.geometry.Triangle`などを用いると，重心や内心，外心，九点円などの三角形が潜在的にもつ様々なデータを自動で計算してくれるので非常に便利です（そこはコンピュータまかせ；）．最終的に，次のような図をプロットするのがゴールです．

![npc.png](https://qiita-image-store.s3.amazonaws.com/0/191414/0644483c-6143-bcae-6ffd-45611488a34d.png "npc.png")



## matplotlibちょい入門

まずは必要なモジュールをインポート．


```python
from sympy.geometry import Point, Circle, Triangle, Segment, Line
import numpy as np
import matplotlib.pyplot as plt
```

### FigureとAxes

matplotlibについての解説はネットにいくらでもありますので，詳しくはそちらをご覧ください（私もそこまで詳しくありません）．今回用いる最低限のものだけ説明します．matplotlibにはFigureとAxesという概念があるのですが，簡単にいうと，`Axes`という座標平面の集まりが`Figure`という理解で大丈夫だと思います．


```python
fig = plt.figure() #Figureオブジェクトを生成
ax1 = fig.add_subplot(1, 2, 1) # figureを1行2列に分けた時の1番目の座標平面を生成
ax2 = fig.add_subplot(1, 2, 2) # figureを1行2列に分けた時の2番目の座標平面を生成
plt.show()
```


![output_5_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/e1e523b9-e25d-1896-c049-553eaba81305.png "output_5_0.png")



`fig.add_subplot`関数によって生成された`Axes`オブジェクトである`ax1`，`ax2`に折れ線や円などを描写していくという流れになります．当面は`Axes`オブジェクトは1つで十分なので，次の記述をベースにしてきます．


```python
fig = plt.figure() #Figureオブジェクトを生成
ax = fig.add_subplot(1, 1, 1) # figureを1行1列に分けた時の1番目の座標平面を生成
plt.show()
```


![output_7_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/ddb83bbb-66f3-6e43-39e8-4db1b27236e5.png "output_7_0.png")



軸の範囲を指定したり，$x$ 軸と $y$ 軸の目盛りを揃えたり，グリッドを追加してみたり，画像として保存したり．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_xlim(-3, 5) #-3<=x<=5 の範囲に設定
ax.set_ylim(-2, 4) #-2<=y<=4 の範囲に設定
ax.set_aspect('equal') #目盛りを揃える
ax.grid() #グリッドを表示
plt.savefig('fig.svg') #figureを画像として保存
plt.show()
```


![output_9_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/d7d10ca4-fc1c-f8fa-0bc8-0b85a12fae49.png "output_9_0.png")



### 折れ線を利用して線分や点を描写しよう

例えば，3点 $(a_1, b_1), (a_2, b_2), (a_3, b_3)$ を順に結んでできる折れ線は次のようにします．

```python
ax.plot([a_1, a_2, a_3], [b_1, b_2, b_3])
```

つまり，$x$ 座標を先に並べて，それに対応する $y$ 座標を順に並べるという引数の取り方をしているのです．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.grid()

ax.plot([1, 2, 3], [3, 2, 5]) #(1, 3), (2, 2), (3, 5)を順に結んでできる折れ線

plt.show()
```


![output_11_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/11fe4a2c-cb43-8893-0aa3-fe3441674f99.png "output_11_0.png")



ここで

```python
p1 = (1, 3)
p2 = (2, 2)
p3 = (3, 5)
```
という3点（タプル）に対して p1，p2，p3 を結んでできる折れ線をプロットするのに`plot(p1, p2, p3)`としてはValue Errorとなります．このような記法を実現するには，[`zip`関数](https://docs.python.jp/3/library/functions.html#zip)の出番です．`zip`関数は
$$\begin{align}
&[a_1, a_2,\ldots,a_n]\\\\
&[b_1, b_2,\ldots,b_n]
\end{align}$$
という複数のリスト（あるいはタプルなど）に対して
$$
[(a_1, b_1), (a_2, b_2),\ldots,(a_n, b_n)]
$$
という1つのリストを生成します．以下は`zip`関数の使用例です．


```python
list(zip([1, 2, 3], ['a', 'b', 'c']))
```




    [(1, 'a'), (2, 'b'), (3, 'c')]




```python
list(zip([1, 2, 3], ['a', 'b', 'c'], ['A', 'B', 'C']))
```




    [(1, 'a', 'A'), (2, 'b', 'B'), (3, 'c', 'C')]




```python
p1 = (1, 3)
p2 = (2, 2)
p3 = (3, 5)
list(zip(p1, p2, p3)) # plot関数用の引数ができる
```




    [(1, 2, 3), (3, 2, 5)]



`zip`関数で次のようなプロットが可能になります．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.grid()

p1 = (1, 3)
p2 = (2, 2)
p3 = (3, 5)

ax.plot(*zip(p1, p2, p3)) #p1, p2, p3を結んでできる折れ線
plt.show()
```


![output_17_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/36bd669a-4020-7a41-9316-34cec0e3f00e.png "output_17_0.png")



引数に`*`を付けていますが，これはリストの中身を引数として扱うための記法です．例えば

```python
x = [x_1, x_2, x_3]
```
のとき，`f(*x)`は`f([x_1, x_2, x_3])`ではなく，

```python
f(x_1, x_2, x_3)
```
の意味になります．


```python
def mysum(a, b, c):
    return a + b + c

#以下はすべて同じ結果
print(mysum(1, 2, 3))
print(mysum(*[1, 2, 3]))
print(mysum(*[1, 2], 3))
```

    6
    6
    6


話を元に戻します．`zip`に突っ込むのは`Point`でもOK．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.grid()

p1 = Point(1, 3)
p2 = Point(2, 2)
p3 = Point(3, 5)
ax.plot(*zip(p1, p2, p3)) #p1, p2, p3を結んでできる折れ線

plt.show()
```


![output_21_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/d97a92d3-cc8f-ade3-2901-278c146ed4f1.png "output_21_0.png")



これならかなり直感的な表記だと思います．始点と終点を一致させることで，*一応*多角形を図示できます．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.grid()

p1 = Point(1, 3)
p2 = Point(2, 2)
p3 = Point(3, 5)
ax.plot(*zip(p1, p2, p3, p1)) #p1, p2, p3を頂点とする三角形

plt.show()
```


![output_23_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/fb1b8a70-c690-62a1-f5fc-57fba12f9bd2.png "output_23_0.png")



線種を破線にするには，引数に`'--'`を追加します．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.grid()

p1 = Point(1, 3)
p2 = Point(2, 2)
p3 = Point(3, 5)
ax.plot(*zip(p1, p2, p3, p1), '--') #線を破線に

plt.show()
```


![output_25_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/a4f74776-acd9-9c10-01c5-8bf6094e9bf0.png "output_25_0.png")



点だけを図示するには，`'o'`などを引数にします．





```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.grid()

p1 = Point(1, 3)
p2 = Point(2, 2)
p3 = Point(3, 5)
ax.plot(*zip(p1, p2, p3, p1), 'o') #頂点のみを表示

plt.show()
```


![output_27_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/1accc365-588b-40c3-47ea-0e0b502820be.png "output_27_0.png")



合わせ技も可能．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.grid()

p1 = Point(1, 3)
p2 = Point(2, 2)
p3 = Point(3, 5)

ax.plot(*zip(p1, p2, p3, p1), 'ro--') #rは赤色の意味

plt.show()
```


![output_29_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/e6408786-399d-d75d-43bd-409edb03cc40.png "output_29_0.png")



線の装飾等などについて詳しくは[matplotlib入門 - プロット方式を変更](http://yubais.net/doc/matplotlib/introduction.html#id4)や[markerの種類](http://matplotlib.org/api/markers_api.html)などをご覧ください．

### 線分や点の描写

2点 $\mathrm{P}_1(1, 3), \mathrm{P}_2(3, 4)$ に対して線分 $\mathrm{P_1P_2}$ を描写するには次のようにすればOK．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.grid()

p1 = Point(1, 3)
p2 = Point(3, 4)
ax.plot(*zip(p1, p2))

plt.show()
```


![output_32_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/16e864b9-8204-12bf-b0af-7754a3e2a404.png "output_32_0.png")



点の描写もplot関数で行えます．点 $(1, 3)$ を描写するのに

```python
p1 = Point(1, 3)
plt.plot(p1[0], p1[1], 'o')
```
とするよりは，やはり`zip`関数を使うとスッキリします（これが果たして可読性が高いといえるのかどうかは分かりませんが……）．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.grid()

p1 = Point(1, 3)
ax.plot(*zip(p1), 'o')

plt.show()
```


![output_34_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/a5c5b4f6-92f2-e57f-9576-e42137a89250.png "output_34_0.png")



### Segmentオブジェクトから線分を図示しよう


```python
P = Point(0, 0)
Q = Point(2, 1)
PQ = Segment(P, Q)
```

で，$\mathrm{P}(0, 0), \mathrm{Q}(2, 1)$ を端点とする線分 $\mathrm{PQ}$ を`Segment`オブジェクトとして考えることができます．`Segment`オブジェクトは`.args`に端点が入ってます．


```python
PQ.args
```




    (Point2D(0, 0), Point2D(2, 1))



したがって `ax.plot(*zip(*PQ.args)`で線分 $\mathrm{PQ}$ を図示できます．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.grid()

ax.plot(*zip(*PQ.args))

plt.show()
```


![output_40_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/22b272a9-a64f-1751-2b42-914d851563dc.png "output_40_0.png")



`ax.plot(*zip(P, Q)`でも同じものが図示できますが，`Segment`オブジェクトから直でプロットできるこの方法も便利です．

## 三角形を図示しよう

`sympy.geometry`の`Triangle`で三角形を定義します．引数は3点を与える方法と辺または角を与える方法があります．また面の概念をもつ図形は`Patch`というクラスに位置付けられており，折れ線とは異なり`add_patch`関数で描写します．三角形の描写には，多角形を表す`Polygon`を使います．引数に頂点の集まりをタプルなどでセットします．

### 3点を指定する方法


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.grid()

p1 = Point(1, 3)
p2 = Point(2, 2)
p3 = Point(3, 5)
t = Triangle(p1, p2, p3)

ax.add_patch(plt.Polygon(t.vertices)) #t.vertices は t の頂点のタプル

plt.show()
```


![output_44_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/06358a33-b142-ff6b-566a-c7d70f9e2ba5.png "output_44_0.png")



何も表示されませんでした．`plot`関数でプロットすると，$x$ 軸や $y$ 軸の値の範囲を自動で設定してくれますが，`add_patch`関数で図形を描写した場合は，デフォルトのものが使われるようです．




```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_xlim(0, 4) #xの範囲を設定
ax.set_ylim(1, 6) #yの範囲を設定
ax.set_aspect('equal')
ax.grid()

p1 = Point(1, 3)
p2 = Point(2, 2)
p3 = Point(3, 5)
t = Triangle(p1, p2, p3)

ax.add_patch(plt.Polygon(t.vertices)) #t.vertices は t の頂点のタプル

plt.show()
```


![output_46_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/89a3cd85-3f6a-78e7-06c2-5730d8bc634e.png "output_46_0.png")



頂点もプロットしておけば一々 $x$ 軸や $y$ 軸の範囲を指定しなくても勝手にやってくれます．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)

ax.set_aspect('equal')
ax.grid()

p1 = Point(1, 3)
p2 = Point(2, 2)
p3 = Point(3, 5)
t = Triangle(p1, p2, p3)

ax.add_patch(plt.Polygon(t.vertices))
ax.plot(*zip(*t.vertices), 'ro') #頂点を赤でプロット

plt.show()
```


![output_48_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/0b0b286a-c056-a5db-bb2c-56eb044a2555.png "output_48_0.png")



### 3辺を指定する方法

`Triangle`には3辺の長さを指定する方法があり，こっちの方が使い勝手が良いです．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.grid()

t = Triangle(sss=(3, 4, 5))

ax.add_patch(plt.Polygon(t.vertices))
ax.plot(*zip(*t.vertices), 'ro')

plt.show()
```


![output_50_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/238944ea-d0f2-746b-6b6a-2ba0c2a01776.png "output_50_0.png")



この場合，原点からスタートして，$x$ 軸の正の方向へ，辺の長さが$3, 4, 5$ となるように三角形が描かれています．

### 2辺とその挟む角を指定する方法

`sas`引数に辺(side)，角(angle)，辺(side)の順でセットします．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.grid()

t = Triangle(sas=(4, 60, 2))

ax.add_patch(plt.Polygon(t.vertices))
ax.plot(*zip(*t.vertices), 'ro')

plt.show()
```


![output_54_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/33bfb71f-12ed-063e-b3c7-f5cc0aa95a3d.png "output_54_0.png")



角度を指定した角が原点にきます．

### 1辺その両端の角を指定する方法

`asa`引数に角(angle)，辺(side)，角(angle)の順でセットします．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.grid()

t = Triangle(asa=(45, 3, 60))

ax.add_patch(plt.Polygon(t.vertices))
ax.plot(*zip(*t.vertices), 'ro')

plt.show()
```


![output_58_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/90911e0c-9ada-d2d8-ba01-e77e9449bff4.png "output_58_0.png")



原点を左端とする $x$ 軸の正の部分に指定した辺がきます．

### 三角形を回転させよう

`sss`引数は便利ですが，必ず底辺が $x$ 軸上に来てしまうので，それが嫌な場合は`rotate`メソッドで回転させましょう．弧度法で指定するので，あらかじめ`pi`をインポートします．


```python
from sympy import pi

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.grid()

t = Triangle(sss=(3, 4, 5)).rotate(pi/2) # 90度回転

ax.add_patch(plt.Polygon(t.vertices))
ax.plot(*zip(*t.vertices), 'ro')

plt.show()
```


![output_61_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/7418105e-680e-9244-9db0-974e84dacb1b.png "output_61_0.png")



### 塗り潰しや辺の色など

塗り潰しをなくすには`fill`引数を`False`にします．


```python
from sympy import pi

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.grid()

t = Triangle(sss=(5, 6, 7))

ax.add_patch(plt.Polygon(t.vertices, fill=False))
ax.plot(*zip(*t.vertices), 'ro')

plt.show()
```


![output_64_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/58e982c6-e752-d326-24e4-e608c7c2d37e.png "output_64_0.png")



線の色や塗り潰しの色の変更，半透明にするには，次のようにします．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.set_axis_off() #軸の非表示

t = Triangle(sss=(5, 6, 7))

ax.add_patch(plt.Polygon(t.vertices,
                         facecolor="orange", #塗り潰しの色
                         linewidth=2, #辺の長さ
                         edgecolor="green", #辺の色
                         alpha=.5 #半透明 0=透明 1=不透明
                        ))
ax.plot(*zip(*t.vertices), 'ro')

plt.show()
```


![output_66_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/c846b691-430d-e6e0-beff-489ff9c6487f.png "output_66_0.png")



### 変数の設定

pythonの命名規則に反するかもしれませんが，数学の伝統にならい，今後，次のような置き方を標準にします．


```python
ABC = Triangle(sss=(5, 6, 7))
A, B, C = ABC.vertices #頂点
AB, BC, CA = Segment(A, B), Segment(B, C), Segment(C, A) #辺の設定 右辺は ABC.sidesと同等
a, b, c = BC.length, CA.length, AB.length #辺の長さ
opposides = { #頂点に対する対辺(opposite side)
    A: BC,
    B: CA,
    C: AB
}
```


```python
print(AB,BC,CA)
print(*ABC.sides)
```

    Segment2D(Point2D(0, 0), Point2D(5, 0)) Segment2D(Point2D(19/5, 12*sqrt(6)/5), Point2D(5, 0)) Segment2D(Point2D(0, 0), Point2D(19/5, 12*sqrt(6)/5))
    Segment2D(Point2D(0, 0), Point2D(5, 0)) Segment2D(Point2D(19/5, 12*sqrt(6)/5), Point2D(5, 0)) Segment2D(Point2D(0, 0), Point2D(19/5, 12*sqrt(6)/5))


`opposide`は造語で<i>opposite side</i>（対辺）の略です．頂点 $\mathrm{A}$ の対辺すなわち  $\mathrm{BC}$ を，`A`をキーとしてアクセスできるととても便利です．

この設定だと，例えば，$\mathrm{AB}$ の中点は `AB.midpoint` で得ることができます．よって中点を結んでできる三角形を次のようにして図示することができます．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.set_axis_off()

ax.add_patch(plt.Polygon((A,B,C), fill=False))
ax.plot(*zip(A,B,C), 'ro')

ax.add_patch(plt.Polygon((BC.midpoint,CA.midpoint,AB.midpoint), fill=False))
ax.plot(*zip(BC.midpoint,CA.midpoint,AB.midpoint), 'ro')

plt.show()
```


![output_71_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/75311e00-ffbb-6c99-fb86-b86d871cea1d.png "output_71_0.png")



### 文字列を表示しよう

三角形の頂点の脇に $\mathrm{A}, \mathrm{B}, \mathrm{C}$ などの文字を置きたくなります．それには，`text`関数を用います．次の例は，点 $(1, 2)$ に $\mathrm{A}$ という文字列を配置します．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.grid()
ax.text(1, 2, '$\mathrm{A}$') #x座標，y座標，内容の順（TeX表記可）
ax.plot(*zip(Point(1, 2)), 'o')
plt.show()
```


![output_73_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/cb72c218-c1f0-bef8-065e-f6e04277dee0.png "output_73_0.png")



点に近いのが少し気になります．`text`には水平方向の位置を表す`ha`引数（`left/center/right`）と，垂直方向の位置を表す`va`引数（`top/center/bottom/baseline`）があります．例えば，点 $(1, 2)$ の左下に $\mathrm{A}$ を配置するには`ha='right' va='top'`とセットします．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.grid()

ax.text(1, 2, '$\mathrm{A}$', ha='right', va='top')
ax.plot(*zip(Point(1, 2)), 'o')
plt.show()
```


![output_75_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/5581ac14-89fa-807a-0bbf-d76d98dd595c.png "output_75_0.png")



直感に反する気もしますが，ようは文字列 $\mathrm{A}$ から見て点 $(1, 2)$ が右上にあるようにセットする，という意味だと思います（この作業は
emathでちまちま点 `\put\A[s]{$1$}` などとやっていたのを思い出させます）．

とりあえず，全パターンに対してどういうプロットがなされるかを次で示しました．やはり点とテキストが近いです．もう少し離す方法はないのでしょうか（`text`の座標をちょっとずらすなどは思いつきますが，座標の目盛りなどに依存しない技が欲しいです）．


```python
from itertools import product

fig = plt.figure(figsize=(8, 8)) #画像サイズを大きめに
ax = fig.add_subplot(1, 1, 1)
ax.set_axis_off()
ax.set_aspect('equal')
ax.set_xlim(-1, 3)
ax.set_ylim(0, 13)

has = ['left', 'center', 'right']
vas = ['top', 'center', 'bottom', 'baseline']

for i, (ha, va) in enumerate(product(has, vas)):
    ax.plot(*zip([0, 12-i]), 'o')
    ax.text(0, 12-i, '$\mathrm{A}$', ha=ha, va=va)
    ax.text(.5, 12-i, "ha='{ha}', va='{va}'".format(ha=ha, va=va), va='center')

plt.savefig('hava.svg') #ブラウザ上で見にくければ，svg画像をご覧ください．
plt.show()
```


![output_77_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/edc9acf5-d74a-3cf2-aa3b-8c614509b5d3.png "output_77_0.png")



さて，本題です．一応，次のようにして三角形に頂点の名前が配置できます．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.set_axis_off()

ax.add_patch(plt.Polygon((A, B, C), fill=False))
ax.plot(*zip(*t.vertices), 'ro')

plt.text(*A, '$\mathrm{A}$', ha='right', va='top')
plt.text(*B, '$\mathrm{B}$', ha='left', va='top')
plt.text(*C, '$\mathrm{C}$', ha='left', va='bottom')

plt.show()
```


![output_79_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/06c85a5b-237d-dfe7-be31-0a30efbecbb7.png "output_79_0.png")



余談ですが，高校の数学では，頂点を表すアルファベットはローマン体で書くことが多いのですが，欧米だと斜体がほとんどのように思えます．頂点 $\mathrm{A}$ の角度，すなわち $\angle\mathrm{A}$を斜体の $A$ で表すので，点としての $\mathrm{A}$ と，数を表す $A$ で区別しているものと思われます．

## 三角形の五心を図示しよう

いよいよコアの部分に突入します．改めて，変数の設定を確認（さっきと同様です）．


```python
ABC = Triangle(sss=(5, 6, 7))
A, B, C = ABC.vertices #頂点
AB, BC, CA = Segment(A, B), Segment(B, C), Segment(C, A) #辺の設定 右辺は ABC.sidesと同等
a, b, c = BC.length, CA.length, AB.length #辺の長さ
opposides = { #頂点に対する対辺(opposite side)
    A: BC,
    B: CA,
    C: AB
}
```

### 重心を図示しよう

重心の定義は，**中線の交点**です．中線とは，頂点とその対辺の中点を結んでできる線分のことです．より具体的にいえば，$\triangle\mathrm{ABC}$ に対して，$\mathrm{A}$ と $\mathrm{BC}$ の中点を結んでできる線分は中線の一つです．`Segment`には`midpoint`属性があります．その名の通り，線分の中点を表します．よって`ax.plot(*zip(A, BC.midpoint))`で中線が引けます．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.set_axis_off()

ax.add_patch(plt.Polygon((A, B, C), fill=False))

G = ABC.centroid #重心を取得
ax.plot(*zip(A, BC.midpoint)) #AとBCの中点を結ぶ
ax.plot(*zip(B, CA.midpoint)) #BとCAの中点を結ぶ
ax.plot(*zip(C, AB.midpoint)) #CとABの中点を結ぶ
ax.plot(*zip(G), 'o') #重心をプロット

plt.text(*A, '$\mathrm{A}$', ha='right', va='top')
plt.text(*B, '$\mathrm{B}$', ha='left', va='top')
plt.text(*C, '$\mathrm{C}$', ha='left', va='bottom')
plt.text(*G, '$\mathrm{G}$', ha='left', va='top')

plt.show()
```


![output_84_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/816466d9-56d3-fbc3-7c0b-4104602789e4.png "output_84_0.png")



対辺を表す`opposides`を用いると`for`ループで中線が引けるので，あとで一括で太さなどを変えることが簡単にできます．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.set_axis_off()

ax.add_patch(plt.Polygon((A, B, C), fill=False))

G = ABC.centroid #重心を取得
for X in [A, B, C]: #中線のプロット
    ax.plot(*zip(X, opposides[X].midpoint), linestyle='--', color='blue', linewidth=1)
ax.plot(*zip(G), 'o') #重心をプロット

ax.text(*A, '$\mathrm{A}$', ha='right', va='top')
ax.text(*B, '$\mathrm{B}$', ha='left', va='top')
ax.text(*C, '$\mathrm{C}$', ha='left', va='bottom')
ax.text(*G, '$\mathrm{G}$', ha='left', va='top')

plt.show()
```


![output_86_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/e0f8c26d-20d5-2dae-2035-5392feb1682d.png "output_86_0.png")



#### 重心の公式

ところで，重心をプロットするのに既に出来上がっているプロパティを用いましたが，高校でも習う通り，$\triangle\mathrm{ABC}$ の重心を $\mathrm{G}$ とすると
$$
\overrightarrow{\mathrm{OG}}=\frac{\overrightarrow{\mathrm{OA}}+\overrightarrow{\mathrm{OB}}+\overrightarrow{\mathrm{OC}}}{3}
$$
で表せます．また，`Point`オブジェクトはベクトルのように加法や実数倍が可能ですので，`G=(A+B+C)/3`で重心が簡単に求まります．


```python
ABC=Triangle(sss=(5,6,7))
A, B, C = ABC.vertices
G = (A+B+C)/3 #重心を計算
G == ABC.centroid
```




    True



### 垂心を図示しよう

垂心の定義は，**頂点から対辺（またはその延長）へ引いた垂線の交点**です．この垂線のことを英語では<i>altitude</i>というそうです（辞書には「高さ」とでます）．`Triangle`オブジェクトには既にこれら垂線が頂点をキーとする辞書として用意されています．


```python
alts = ABC.altitudes #垂線の取得
alts
```




    {Point2D(0, 0): Segment2D(Point2D(0, 0), Point2D(24/5, 2*sqrt(6)/5)),
     Point2D(5, 0): Segment2D(Point2D(361/245, 228*sqrt(6)/245), Point2D(5, 0)),
     Point2D(19/5, 12*sqrt(6)/5): Segment2D(Point2D(19/5, 0), Point2D(19/5, 12*sqrt(6)/5))}



`alts[A]`は頂点 $\mathrm{A}$ からの垂線を表す線分が返ってくるので，[さっき説明した](#Segmentオブジェクトから線分を図示しよう)`ax.plot(*zip(*alts[A].args))`で垂線が図示できます．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.set_axis_off()
ax.add_patch(plt.Polygon((A, B, C), fill=False))

H=ABC.orthocenter #垂心の取得
alts = ABC.altitudes #垂線の取得
for X in [A, B, C]: #垂線のプロット
    ax.plot(*zip(*alts[X].args))
ax.plot(*zip(H), 'o') #垂心をプロット

plt.text(*A, '$\mathrm{A}$', ha='right', va='top')
plt.text(*B, '$\mathrm{B}$', ha='left', va='top')
plt.text(*C, '$\mathrm{C}$', ha='left', va='bottom')
plt.text(*H, '$\mathrm{H}$', ha='left', va='bottom')

plt.show()
```


![output_93_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/92b0daa6-976c-33be-2272-222c0873ce5a.png "output_93_0.png")



### 線分ではなく直線を図示しよう

外心の図示に進む前に，直線のプロットを考えます．数学的に直線とは始点と終点が存在しない無限に伸びた線のことです．そもそもmatplotlibはヒストグラムや散布図など，有限のデータを扱うためのものだと思っているので，無限に長い直線なんてのは本来の用途とは外れる気がするのですが，今回は擬似的に，**とても長い線分**として直線をプロットすることにします．

#### 直線のベクトル方程式

今したいのは，異なる2点 $\mathrm{P}, \mathrm{Q}$ が与えられたときに，線分 $\mathrm{PQ}$ を延長したとても長い線分をプロットすることです．高校で習った通り，直線 $\mathrm{PQ}$ のベクトル方程式は
$$
\overrightarrow{\mathrm{OX}} = \overrightarrow{\mathrm{OP}} + t\overrightarrow{\mathrm{PQ}} \quad(-\infty < t <\infty)
$$
です．$t$ の範囲をいじれば長さを変えることができます．例えば線分 $\mathrm{PQ}$ のベクトル方程式は
$$
\overrightarrow{\mathrm{OX}} = \overrightarrow{\mathrm{OP}} + t\overrightarrow{\mathrm{PQ}} \quad(0 \leqq t \leqq 1)
$$
です．線分 $\mathrm{PQ}$ を $\mathrm{Q}$ の方向へ$2$倍にした線分は
$$
\overrightarrow{\mathrm{OX}} = \overrightarrow{\mathrm{OP}} + t\overrightarrow{\mathrm{PQ}} \quad(0 \leqq t \leqq 2)
$$
になります．

一般に $t_m<t_M$ なる実数に対して
$$
\overrightarrow{\mathrm{OX}} = \overrightarrow{\mathrm{OP}} + t\overrightarrow{\mathrm{PQ}} \quad(t_m \leqq t \leqq t_M)
$$
と表せる線分について考えましょう．
$$
\overrightarrow{\mathrm{PQ}}=\overrightarrow{\mathrm{OQ}}-\overrightarrow{\mathrm{OP}}
$$
ですから，この線分の始点の位置ベクトルは $t=t_m$ を代入した
$$
\overrightarrow{\mathrm{OP}} + t_m(\overrightarrow{\mathrm{OQ}}-\overrightarrow{\mathrm{OP}})
$$
で，終点の位置ベクトルは $t=t_M$ を代入した
$$
\overrightarrow{\mathrm{OP}} + t_M(\overrightarrow{\mathrm{OQ}}-\overrightarrow{\mathrm{OP}})
$$
で表せます．これを python で表現すると，始点は`p+(q-p)*t_m`，終点は`p+(q-p)*t_M`となります．実際に
```python
p = Point(0, 0)
q = Point(2, 1)
```
という設定で $\mathrm{PQ}$ を延長した線分をプロットしてみましょう．次は
$$
\overrightarrow{\mathrm{OX}} = \overrightarrow{\mathrm{OP}} + t\overrightarrow{\mathrm{PQ}} \quad(0 \leqq t \leqq 2)
$$
のプロット例です．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')

p = Point(0, 0)
q = Point(2, 1)

t_m = 0
t_M = 2

start = p+(q-p)*t_m #始点
end = p+(q-p)*t_M #終点

ax.plot(*zip(p, q), linewidth=4) #線分PQは太く
ax.plot(*zip(start, end), linewidth=1)

ax.text(*p, '$\mathrm{P}$', ha='left', va='top')
ax.text(*q, '$\mathrm{Q}$', ha='left', va='top')

plt.show()
```


![output_96_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/59fe4db4-3dfd-c61e-3642-ae16882cb75d.png "output_96_0.png")



確かに $\mathrm{PQ}$ を2倍にした線分がプロットされました．次は
$$
\overrightarrow{\mathrm{OX}} = \overrightarrow{\mathrm{OP}} + t\overrightarrow{\mathrm{PQ}} \quad(-2 \leqq t \leqq 3)
$$
を図示してみます．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')

p = Point(0, 0)
q = Point(2, 1)

t_m = -2
t_M = 3

start = p+(q-p)*t_m #始点
end = p+(q-p)*t_M #終点

ax.plot(*zip(p, q), linewidth=4) #線分PQは太く
ax.plot(*zip(start, end), linewidth=1)

ax.text(*p, '$\mathrm{P}$', ha='left', va='top')
ax.text(*q, '$\mathrm{Q}$', ha='left', va='top')

plt.show()
```


![output_98_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/f77688af-3214-d6d2-479d-d177de8a5731.png "output_98_0.png")



程度がよく分かりませんが，`t_m=-100`，`t_M=100`とすればほぼ無限に長い直線といえそうです．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')

p = Point(0, 0)
q = Point(2, 1)

t_m = -100
t_M = 100

start = p+(q-p)*t_m #始点
end = p+(q-p)*t_M #終点

ax.plot(*zip(p, q), linewidth=4) #線分PQは太く
ax.plot(*zip(start, end), linewidth=1)

ax.text(*p, '$\mathrm{P}$', ha='left', va='top')
ax.text(*q, '$\mathrm{Q}$', ha='left', va='top')

plt.show()
```


![output_100_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/60e279d0-ecd7-fd87-63b7-9521f4c33246.png "output_100_0.png")



実際に使うときは，軸の範囲を狭めます．


```python
plt.show()
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.set_xlim(-10, 10) #xの範囲を絞る
ax.set_ylim(-10, 10) #yの範囲を絞る

p = Point(0, 0)
q = Point(2, 1)

t_m = -100
t_M = 100

start = p+(q-p)*t_m #始点
end = p+(q-p)*t_M #終点

ax.plot(*zip(p, q), linewidth=4) #線分PQは太く
ax.plot(*zip(start, end), linewidth=1)

ax.text(*p, '$\mathrm{P}$', ha='left', va='top')
ax.text(*q, '$\mathrm{Q}$', ha='left', va='top')

plt.show()
```


![output_102_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/1f853c45-e95e-86f0-a87c-f4fa763fdc6d.png "output_102_0.png")



関数にして扱いやすくしましょう．


```python
def side_of_line(p, q, mint=-100, maxt=100):
    return zip(p+(q-p)*mint, p+(q-p)*maxt)

def side_of_ray(p, q, maxt=100):
    return zip(p, p+(q-p)*maxt)
```

sideは端点，という意味のつもりで使ってます．ついでに半直線（片方だけ無限に長い線．英語で<i>ray</i>）も定義しておきました．次の例は，$\mathrm{P}(0, 0), \mathrm{Q}(2, 1), \mathrm{R}(0, 1), \mathrm{S}(2, 2)$ に対して直線 $\mathrm{PQ}$ と，半直線 $\mathrm{RS}$ をプロットする例です．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)

p = Point(0, 0)
q = Point(2, 1)
r = Point(0, 1)
s = Point(2, 2)

ax.plot(*side_of_line(p, q))
ax.plot(*side_of_ray(r, s))

ax.text(*p, '$\mathrm{P}$', ha='left', va='top')
ax.text(*q, '$\mathrm{Q}$', ha='left', va='top')
ax.text(*r, '$\mathrm{R}$', ha='right', va='bottom')
ax.text(*s, '$\mathrm{S}$', ha='right', va='bottom')

ax.plot(*zip(p, q, r, s), '.')

plt.show()
```


![output_106_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/30209d8b-4013-5eef-25a3-4064fad0cbdb.png "output_106_0.png")



### 外心と外接円を図示しよう

外心は*外接円の中心*といっても同じことですが，伝統的には**各辺の垂直二等分線の交点**です．`Triangle`オブジェクトには当然のことながら垂直二等分線も用意されています．線分 `AB` の垂直二等分線は`AB.perpendicular_bisector()`で得られます．これをSegmentオブジェクトと同様にしてプロットをすると一つ問題が発生します．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.set_axis_off()

ax.add_patch(plt.Polygon((A, B, C), fill=False))

ax.plot(*zip(*AB.perpendicular_bisector().args)) #ABの垂直二等分線をプロット？
ax.plot(*zip(*BC.perpendicular_bisector().args)) #BCの垂直二等分線をプロット？
ax.plot(*zip(*CA.perpendicular_bisector().args)) #CAの垂直二等分線をプロット？
ax.plot(*ABC.circumcenter, 'o') #外心をプロット

ax.text(*A, '$\mathrm{A}$', ha='right', va='top')
ax.text(*B, '$\mathrm{B}$', ha='left', va='top')
ax.text(*C, '$\mathrm{C}$', ha='left', va='bottom')

plt.show()
```


![output_108_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/a5583a26-2e34-209a-5305-e8c30e7b05af.png "output_108_0.png")



`AB.perpendicular_bisector().args`によって得られる直線上の2点がどういう風に設定されているのかは不明（片方は$\mathrm{AB}$の中点）ですが，やはり短すぎです．そこで，`side_of_line`関数で200倍にしましょう．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.set_axis_off()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)

ax.add_patch(plt.Polygon((A, B, C), fill=False))

ax.plot(*side_of_line(*AB.perpendicular_bisector().args)) #ABの垂直二等分線をプロット
ax.plot(*side_of_line(*BC.perpendicular_bisector().args)) #BCの垂直二等分線をプロット
ax.plot(*side_of_line(*CA.perpendicular_bisector().args)) #CAの垂直二等分線をプロット

ax.plot(*ABC.circumcenter, 'o') #外心

ax.text(*A, '$\mathrm{A}$', ha='right', va='top')
ax.text(*B, '$\mathrm{B}$', ha='left', va='top')
ax.text(*C, '$\mathrm{C}$', ha='left', va='bottom')

plt.show()
```


![output_110_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/cadaf89c-71ef-5b33-c82b-04346325b57a.png "output_110_0.png")



うまくいきました．もう少しスリム化して，外接円などを足したものが次です．


```python
ax = plt.gca()
ax.set_xlim((-3,8))
ax.set_ylim((-3,8))
ax.set_axis_off()
ax.set_aspect('equal')

ax.add_patch(plt.Polygon((A, B, C), fill=False))

O=ABC.circumcenter #外心の取得

for X in [A, B, C]:
    ax.plot(*side_of_line(*opposides[X].perpendicular_bisector().args)) #Xの対辺の垂直二等分線をプロット

ax.plot(*O, 'o') #外心をプロット
ax.add_patch(plt.Circle(O, ABC.circumradius, fill=False)) #外接円をプロット

ax.text(*A, '$\mathrm{A}$', ha='right', va='top')
ax.text(*B, '$\mathrm{B}$', ha='left', va='top')
ax.text(*C, '$\mathrm{C}$', ha='left', va='bottom')
ax.text(*O, '$\mathrm{O}$', ha='right', va='top')

plt.show()
```


![output_112_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/c5de0a50-314d-d9fc-07db-cdaaae9ae10e.png "output_112_0.png")



### 内心と内接円を図示しよう

内心は**内角の二等分線の交点**であり，*内接円の中心*でもあります．`ABC.bisectors`に頂点をキーとする角の二等分線を表す`Segment`が入っています．


```python
ax = plt.gca()
ax.set_axis_off()
ax.set_aspect('equal')

ax.add_patch(plt.Polygon((A, B, C), fill=False))

bisectors = ABC.bisectors() #内角の二等分線を取得

I = ABC.incenter #内心の取得
for X in [A, B, C]:
    ax.plot(*zip(*bisectors[X].args)) #角Xの二等分線をプロット
ax.plot(*zip(I), 'o') #内心をプロット
ax.add_patch(plt.Circle(I, ABC.inradius, fill=False)) #内接円をプロット

ax.text(*A, '$\mathrm{A}$', ha='right', va='top')
ax.text(*B, '$\mathrm{B}$', ha='left', va='top')
ax.text(*C, '$\mathrm{C}$', ha='left', va='bottom')
ax.text(*I, '$\mathrm{I}$', ha='left', va='top')

plt.show()
```


![output_114_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/0467a42f-bc2f-53e9-27ed-2f18bf2e9250.png "output_114_0.png")



### 傍心と傍接円を図示しよう

傍心は，外角の二等分線交点として定義される点で少し難しいのですが，ある意味で内心を拡張したような存在です．残念ながら`Triangle`にはそれらしい属性が存在しませんでした．なので自力で作るしかありません．とりあえず，次のコードは無視して出来上がった図を見てください．


```python
ax = plt.gca()
ax.set_axis_off()
ax.set_aspect('equal')
ax.set_xlim((-6, 8))
ax.set_ylim((-5, 8))

ax.add_patch(plt.Polygon((A, B, C), fill=False))

ax.plot(*side_of_line(A, B), color='black') #ABの延長

D = A*2-B

CAD = Triangle(C, A, D) #三角形CAD

ax.plot(*side_of_line(*CAD.bisectors()[A].args), linewidth=2, color='red') #CADのAの二等分線

ax.text(*A, '$\mathrm{A}$', ha='right', va='top')
ax.text(*B, '$\mathrm{B}$', ha='left', va='top')
ax.text(*C, '$\mathrm{C}$', ha='left', va='bottom')
ax.text(*D, '$\mathrm{D}$', ha='right', va='top')

ax.plot(*zip(A,B,C,D),'o')

plt.show()
```


![output_117_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/2de7c729-a88a-6d25-cac4-1b12b5cdf2c7.png "output_117_0.png")



図の $\angle\mathrm{CAD}$ を $\mathrm{A}$ の外角といいます．外角の二等分線とは，図の赤色の直線です．$\mathrm{AB}=\mathrm{AD}$ となるように点 $\mathrm{D}$ をとっているので，
$$
\overrightarrow{\mathrm{OD}} = \overrightarrow{\mathrm{OA}}-\overrightarrow{\mathrm{AB}}
=\overrightarrow{\mathrm{OA}}-(\overrightarrow{\mathrm{OB}}-\overrightarrow{\mathrm{OA}})
=\overrightarrow{\mathrm{2OA}}-\overrightarrow{\mathrm{OB}}
$$
ですから

```python
D = A*2-B
```
で点  $\mathrm{D}$ が作れます．そして，三角形 $\triangle\mathrm{CAD}$ の角 $\angle\mathrm{A}$ の二等分線として外角の二等分線をプロットしているのです．

全ての頂点の外角の二等分線をプロットしてみます．


```python
ax = plt.gca()
ax.set_axis_off()
ax.set_aspect('equal')
ax.set_xlim((-6, 10))
ax.set_ylim((-5, 9))

ax.add_patch(plt.Polygon((A, B, C), fill=False))

ax.plot(*side_of_line(A, B), color='black') #ABの延長
ax.plot(*side_of_line(B, C), color='black') #BCの延長
ax.plot(*side_of_line(C, A), color='black') #CAの延長

ax.plot(*side_of_line(*Triangle(C, A, A*2-B).bisectors()[A].args), linewidth=2, color="red")
ax.plot(*side_of_line(*Triangle(A, B, B*2-C).bisectors()[B].args), linewidth=2, color="red")
ax.plot(*side_of_line(*Triangle(B, C, C*2-A).bisectors()[C].args), linewidth=2, color="red")

ax.text(*A, '$\mathrm{A}$', ha='right', va='top')
ax.text(*B, '$\mathrm{B}$', ha='left', va='top')
ax.text(*C, '$\mathrm{C}$', ha='left', va='bottom')

ax.plot(*zip(A,B,C),'o')

plt.show()
```


![output_119_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/0c6f5aee-e8ba-389c-f528-0686476dea2d.png "output_119_0.png")



これら外角の二等分線の3つの交点が傍心です．sympyは二直線に対して交点を求める関数が用意されていますが，ここは天下り的に傍心の座標を計算します．
三角形 $\triangle\mathrm{ABC}$ を固定するとき，任意の平面上の点 $\mathrm{P}$ は $3$ つの実数 $s, t, u$ を用いて
$$
\overrightarrow{\mathrm{OP}}=\frac{s\overrightarrow{\mathrm{OA}}+t\overrightarrow{\mathrm{OB}}+u\overrightarrow{\mathrm{OC}}}{s+t+u}
$$
という形で表せることが知られています．この3つの実数の組 $(s, t, u)$ を $\mathrm{P}$ の**重心座標**といいます．重心 $\mathrm{G}$ の重心座標は $(1, 1, 1)$ です．実は，内心 $\mathrm{I}$ は
$$
\overrightarrow{\mathrm{OI}}=\frac{a\overrightarrow{\mathrm{OA}}+b\overrightarrow{\mathrm{OB}}+c\overrightarrow{\mathrm{OC}}}{a+b+c}
$$
で表せることが知られているので，重心座標は $(a, b, c)$ です．実は，3つの傍心の重心座標は
$$
(-a, b, c),\ (a, -b, c),\ (a, b, -c)
$$
であることが知られています．よって傍心は
$$
\frac{-a\overrightarrow{\mathrm{OA}}+b\overrightarrow{\mathrm{OB}}+c\overrightarrow{\mathrm{OC}}}{-a+b+c},\ 
\frac{a\overrightarrow{\mathrm{OA}}-b\overrightarrow{\mathrm{OB}}+c\overrightarrow{\mathrm{OC}}}{a-b+c},\ 
\frac{a\overrightarrow{\mathrm{OA}}+b\overrightarrow{\mathrm{OB}}-c\overrightarrow{\mathrm{OC}}}{a+b-c}
$$
で計算できます．そして，傍心を中心とし各辺の延長に接する円，傍接円というものがあり，それらの半径は
$$
\frac{2S}{-a+b+c},\ \frac{2S}{a-b+c},\ \frac{2S}{a+b-c}
$$
で与えられます．ここで $S$ は $\triangle\mathrm{ABC}$ の面積です．これらを総括してプロットしたものの例が次です．


```python
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(1, 1, 1)
ax.set_axis_off()
ax.set_aspect('equal')
ax.set_xlim((-15, 15))
ax.set_ylim((-10, 18))

ax.add_patch(plt.Polygon((A, B, C), fill=False))

for seg in [AB, BC, CA]: #線分の延長
    ax.plot(*side_of_line(*seg.args))

ax.plot(*side_of_line(*Triangle(B, A, A*2-C).bisectors()[A].args), linewidth=.75)
ax.plot(*side_of_line(*Triangle(C, B, B*2-A).bisectors()[B].args), linewidth=.75) #外角の二等分線
ax.plot(*side_of_line(*Triangle(A, C, C*2-B).bisectors()[C].args), linewidth=.75)

excenters = { #傍心
    A: (-A*a+B*b+C*c)/(-a+b+c),
    B: (A*a-B*b+C*c)/(a-b+c),
    C: (A*a+B*b-C*c)/(a+b-c)
}

exradiuses = { #傍接円の半径 ABC.areaはABCの面積
    A: 2*ABC.area/(-a+b+c),
    B: 2*ABC.area/(a-b+c),
    C: 2*ABC.area/(a+b-c)
}

ax.plot(*zip(*excenters.values()), 'o') #傍心のプロット

for X in [A, B, C]: #傍接円のプロット
    ax.add_patch(plt.Circle(excenters[X], exradiuses[X], fill=False))

ax.text(*A, '$\mathrm{A}$', ha='right', va='top')
ax.text(*B, '$\mathrm{B}$', ha='left', va='top')
ax.text(*C, '$\mathrm{C}$', ha='center', va='bottom')
ax.text(*excenters[A], '$\mathrm{E_A}$', ha='center', va='bottom')
ax.text(*excenters[B], '$\mathrm{E_B}$', ha='left', va='bottom')
ax.text(*excenters[C], '$\mathrm{E_C}$', ha='right', va='center')
    
plt.show()
```


![output_121_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/65626fbe-0e3d-af48-995a-89605ef3b41d.png "output_121_0.png")



### オイラー線

重心と外心と垂心は同一直線上に乗っかっていることが知られています．この直線をオイラー線といいます．確認してみましょう．


```python
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(1, 1, 1)
ax.set_axis_off()
ax.set_aspect('equal')
ax.set_xlim((-3,8))
ax.set_ylim((-3,8))

ax.add_patch(plt.Polygon((A, B, C), fill=False))

G=ABC.centroid #重心の取得
O=ABC.circumcenter #外心の取得
H=ABC.orthocenter #垂心の取得

ax.plot(*side_of_line(O, H)) #オイラー線のプロット
ax.plot(*zip(G, O, H), 'o') #3点のプロット

ax.text(*A, '$\mathrm{A}$', ha='right', va='top')
ax.text(*B, '$\mathrm{B}$', ha='left', va='top')
ax.text(*C, '$\mathrm{C}$', ha='left', va='bottom')
ax.text(*G, '$\mathrm{G}$', ha='right', va='top')
ax.text(*O, '$\mathrm{O}$', ha='right', va='top')
ax.text(*H, '$\mathrm{H}$', ha='right', va='top')

plt.show()
```


![output_123_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/e251293b-487a-4272-8fc8-41a7e673a490.png "output_123_0.png")



ちなみに，$\mathrm{OG}: \mathrm{GH} = 1:2$ であることが知られています．
$$
\overrightarrow{\mathrm{OH}} = \overrightarrow{\mathrm{OA}}+\overrightarrow{\mathrm{OB}}+\overrightarrow{\mathrm{OC}}
$$
であることを示せば重心の公式より
$$
\overrightarrow{\mathrm{OG}} = \frac{\overrightarrow{\mathrm{OA}}+\overrightarrow{\mathrm{OB}}+\overrightarrow{\mathrm{OC}}}{3}
$$
であることより直ちに証明できます．

## 九点円を図示しよう

いよいよ九点円に入ります．

### どの九点？

三角形 $\triangle\mathrm{ABC}$ において，

+ 各辺の中点
+ 頂点から対辺に引いた垂線の足
+ 垂心と頂点の中点

計9つをプロットしましょう．


```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.set_axis_off()

H = ABC.orthocenter #垂心を取得
alts = ABC.altitudes #垂線を取得

nine_points = []

for X in [A, B, C]:
    nine_points.append(opposides[X].midpoint) #対辺の中点
    if X == alts[X].p1: #線分alts[X] の始点alts[X].p1, 終点alts[X].p2のうち，頂点でない方が垂線の足
        nine_points.append(alts[X].p2)
    else:
        nine_points.append(alts[X].p1)
    nine_points.append(Segment(X, H).midpoint) #垂心と頂点の中点

ax.add_patch(plt.Polygon((A, B, C), fill=False))

for X in [A, B, C]:
    ax.plot(*zip(*alts[X].args))
ax.plot(*zip(*nine_points), 'ro')

plt.show()
```


![output_126_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/715a3c56-00bd-b20c-5aab-82cde8307e62.png "output_126_0.png")



確かにこの九つの点は同一の円上にありそうです．なんと九点円は初めから`Triagnle`オブジェクトに入っています．


```python
npc = ABC.nine_point_circle #ABCの九点円
print(npc.center, npc.radius) #九点円の中心と半径
```

    Point2D(63/20, 163*sqrt(6)/240) 35*sqrt(6)/48



```python
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.set_axis_off()

H = ABC.orthocenter #垂心を取得
alts = ABC.altitudes #垂線を取得

nine_points = []

for X in [A, B, C]:
    nine_points.append(opposides[X].midpoint) #対辺の中点
    if X == alts[X].p1: #線分alts[X] の始点alts[X].p1, 終点alts[X].p2のうち，頂点でない方が垂線の足
        nine_points.append(alts[X].p2)
    else:
        nine_points.append(alts[X].p1)
    nine_points.append(Segment(X, H).midpoint) #垂心と頂点の中点

ax.add_patch(plt.Polygon((A, B, C), fill=False))

for X in [A, B, C]:
    ax.plot(*zip(*alts[X].args))
ax.plot(*zip(*nine_points), 'ro')

npc = ABC.nine_point_circle #九点円を取得
ax.add_patch(plt.Circle(npc.center, npc.radius, fill=False, linewidth=1)) #九点円をプロット
ax.plot(*zip(npc.center), 'o', color = 'purple') #九点円の中心をプロット

plt.show()
```


![output_129_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/aa2559bc-4bd9-5e29-c4e3-40c72b886a94.png "output_129_0.png")



### 九点円の中心とオイラー線

九点円の中心は外心と垂心の中点であることが知られています．ということは，九点円の中心と重心，外心，垂心の4点はオイラー上にあることが分かります．


```python
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(1, 1, 1)
ax.set_xlim((-2,7))
ax.set_ylim((-2,7))
ax.set_aspect('equal')
ax.set_axis_off()

ax.add_patch(plt.Polygon((A, B, C), fill=False))

npc = ABC.nine_point_circle #九点円を取得
ax.add_patch(plt.Circle(npc.center, npc.radius, fill=False, linewidth=1)) #九点円をプロット
ax.plot(*zip(npc.center), 'o', color = 'purple') #九点円の中心をプロット

G = ABC.centroid #重心の取得
O = ABC.circumcenter #外心の取得
H = ABC.orthocenter #垂心の取得

ax.plot(*side_of_line(O, H)) #オイラー線のプロット
ax.plot(*zip(G, O, H), 'o') #3点のプロット

ax.text(*A, '$\mathrm{A}$', ha='right', va='top')
ax.text(*B, '$\mathrm{B}$', ha='left', va='top')
ax.text(*C, '$\mathrm{C}$', ha='center', va='bottom')
ax.text(*G, '$\mathrm{G}$', ha='right', va='top')
ax.text(*O, '$\mathrm{O}$', ha='right', va='top')
ax.text(*H, '$\mathrm{H}$', ha='right', va='top')
ax.text(*npc.center, '$\mathrm{N}$', ha='right', va='top')

plt.show()
```


![output_131_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/e6bb9004-904f-1e1d-0444-dd58de9eed06.png "output_131_0.png")



### フォイエルバッハの定理

九点円は内接円に内接し，傍接円に外接しています．


```python
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.set_axis_off()
ax.set_xlim((-1, 8))
ax.set_ylim((-3, 6))

ax.add_patch(plt.Polygon((A, B, C), fill=False))
for seg in [AB, BC, CA]: #線分の延長
    ax.plot(*side_of_line(*seg.args), color='black')

npc = ABC.nine_point_circle #ABCの九点円

I = ABC.incenter

excenters = { #傍心
    A: (-A*a+B*b+C*c)/(-a+b+c),
    B: (A*a-B*b+C*c)/(a-b+c),
    C: (A*a+B*b-C*c)/(a+b-c)
}

exradiuses = { #傍接円の半径 ABC.areaはABCの面積
    A: 2*ABC.area/(-a+b+c),
    B: 2*ABC.area/(a-b+c),
    C: 2*ABC.area/(a+b-c)
}

ax.add_patch(plt.Circle(I, ABC.inradius, fill=False, linewidth=1, color='green')) #内接円のプロット

for X in [A, B, C]: #傍接円のプロット
    ax.add_patch(plt.Circle(excenters[X], exradiuses[X], fill=False, linewidth=1, color='green'))

ax.add_patch(plt.Circle(npc.center, npc.radius, fill=False, linewidth=2, color="red"))

plt.show()
```


![output_133_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/3d6ae160-4339-8328-a296-193e342bdb4b.png "output_133_0.png")



確かに九点円（赤）は，内接円と傍接円と接しています．

まとめると九点円は

+ 中心が重心，外心，垂心と直線で結ばれている
+ 内心，傍心を中心とする内接円，傍接円に接している．

ということが分かり，三角形の五心である，重心，外心，垂心，内心，傍心と関係を持っていることが分かりました．

### 最後に全部詰め込もう


```python
from sympy.geometry import Point, Circle, Triangle, Segment, Line
import numpy as np
import matplotlib.pyplot as plt

def side_of_line(p, q, mint=-100, maxt=100):
    return zip(p+(q-p)*mint, p+(q-p)*maxt)

ABC = Triangle(sss=(5, 6, 7))
A, B, C = ABC.vertices #頂点
AB, BC, CA = Segment(A, B), Segment(B, C), Segment(C, A) #辺の設定 右辺は ABC.sidesと同等
a, b, c = BC.length, CA.length, AB.length #辺の長さ
opposides = {A: BC, B: CA, C: AB} #頂点に対する対辺(opposite side)

fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.set_axis_off()
ax.set_xlim((-10, 14))
ax.set_ylim((-8, 16))

npc = ABC.nine_point_circle #ABCの九点円を取得
N = npc.center #九点円の中心を取得
G = ABC.centroid #重心の取得
O = ABC.circumcenter #外心の取得
H = ABC.orthocenter #垂心の取得
I = ABC.incenter #内心の取得

excenters = { #傍心を計算
    A: (-A*a+B*b+C*c)/(-a+b+c),
    B: (A*a-B*b+C*c)/(a-b+c),
    C: (A*a+B*b-C*c)/(a+b-c)
}

exradiuses = { #傍接円の半径 ABC.areaはABCの面積
    A: 2*ABC.area/(-a+b+c),
    B: 2*ABC.area/(a-b+c),
    C: 2*ABC.area/(a+b-c)
}

alts = ABC.altitudes #垂線を取得

nine_points = [] #九点
for X in [A, B, C]:
    nine_points.append(opposides[X].midpoint) #対辺の中点
    if X == alts[X].p1: #線分alts[X] の始点alts[X].p1, 終点alts[X].p2のうち，頂点でない方が垂線の足
        nine_points.append(alts[X].p2)
    else:
        nine_points.append(alts[X].p1)
    nine_points.append(Segment(X, H).midpoint) #垂心と頂点の中点
ax.add_patch(plt.Polygon((A, B, C), facecolor='#67BE91', alpha=.3)) #三角形ABCをプロット
for seg in ABC.sides:
    ax.plot(*side_of_line(*seg.args), color='#422500', linewidth=.3) #辺の延長をプロット

for X in [A, B, C]:
    seg = opposides[X] #Xの対辺
    ax.plot(*zip(*alts[X].args), color='#00612E', linewidth=.2) #垂線をプロット
    ax.plot(*side_of_line(*seg.perpendicular_bisector().args), '--', color='gray', linewidth=.3) #垂直二等分線をプロット
    ax.plot(*zip(X, seg.midpoint), color='#00612E', linewidth=.2) #中線をプロット
    ax.plot(*zip(*ABC.bisectors()[X].args), '--', color='gray', linewidth=.3) #角Xの二等分線

ax.plot(*side_of_line(*Triangle(B, A, A*2-C).bisectors()[A].args),'--',  color='gray', linewidth=.3)
ax.plot(*side_of_line(*Triangle(C, B, B*2-A).bisectors()[B].args),'--',  color='gray', linewidth=.3) #外角の二等分線
ax.plot(*side_of_line(*Triangle(A, C, C*2-B).bisectors()[C].args),'--',  color='gray', linewidth=.3)

ax.plot(*side_of_line(O, H), '-', linewidth=5, color='#EC6800', alpha=.4)#オイラー線をプロット

ax.add_patch(plt.Circle(I, ABC.inradius, linestyle='--', edgecolor='#9D221B', fill=False, linewidth=.3, alpha=1)) #内接円をプロット
ax.add_patch(plt.Circle(I, ABC.inradius, facecolor='gray', linewidth=0, alpha=.2)) #内接円の内部を塗り潰す

for X in [A, B, C]: #傍接円をプロット
    ax.add_patch(plt.Circle(excenters[X], exradiuses[X], edgecolor='gray', facecolor='#8BA55D', linewidth=0, alpha=.3))

ax.add_patch(plt.Circle(N, npc.radius, edgecolor='#EC6800', linewidth=.4, fill=False)) #九点円をプロット
ax.add_patch(plt.Circle(N, npc.radius, linewidth=0, facecolor='#64954F',fill=True, alpha=.3)) #九点円の内部を塗りつぶす

ax.add_patch(plt.Circle(O, ABC.circumradius, linestyle='--', edgecolor='gray', fill=False, linewidth=.3)) #外接円をプロット

ax.add_patch(plt.Polygon(tuple(excenters.values()), facecolor='#C3D300', alpha=.1)) #傍心を頂点とする三角形を塗り潰し

ax.plot(*zip(G, O, H, N), '.', color='#EF3706', ms=3) #重心，外心，垂心，九点円の中心をプロット
ax.plot(*zip(A, B, C), '.', color='black', ms=3, alpha=.5) #頂点をプロット
ax.plot(*zip(*excenters.values()), '.', color='#074d59', alpha=.5, ms=3) #傍心のプロット
ax.plot(*zip(ABC.incenter), '.', color='#074d59', alpha=.5, ms=3) #内心のプロット
ax.plot(*zip(*nine_points), '.', color='#5E4500', ms=2, alpha=.5) #九点をプロット

plt.savefig('npc.png', dpi=100)
plt.savefig('npc.svg')
plt.show()
```


![output_136_0.png](https://qiita-image-store.s3.amazonaws.com/0/191414/5178c103-e2e1-a3c3-8898-ce8edda24c2d.png "output_136_0.png")



## Triangleを拡張したクラスを作りたい

`Triangle`に用意されていない傍心の座標などをプロパティとしてもつ`Triangle`を継承したクラスを作りたいと思うのですが，どうも，私の生半可な知識ではうまくいきません．


```python
from sympy.geometry import Triangle

class MyTriangle(Triangle):
    @property
    def nine_point_center(self):
        return self.nine_point_circle.center      
```

プロパティを追加するだけなら，これで大丈夫かと思ったのですが，


```python
a, b, c = (0, 0), (2, 0), (1, 1)
t1 = MyTriangle(a, b, c)
print(t1.nine_point_center)
```

    Point2D(1, 1/2)


はうまくいくのですが，


```python
t2 = MyTriangle(sss=(5, 6, 7))
print(t2.nine_point_center)
```


    ---------------------------------------------------------------------------

    AttributeError                            Traceback (most recent call last)

    <ipython-input-65-1d4a6d686027> in <module>()
          1 t2 = MyTriangle(sss=(5, 6, 7))
    ----> 2 print(t2.nine_point_center)
    

    AttributeError: 'Triangle' object has no attribute 'nine_point_center'


だとうまくいきません．`t2`は`Triangle`オブジェクトとして認識されているようです．`__new__()`なども色々試してみましたがうまくいきませんでした．もし解決策をご存知の方がいらっしゃったらご教示お願いします．

## 参考
+ [matplotlib入門](http://bicycle1885.hatenablog.com/entry/2014/02/14/023734)