circumcircle, incircle and excircle of a triangle
概要
三点 \(\mathbf{p_1}, \mathbf{p_2}, \mathbf{p_3}\) から成る三角形の外接円と内接円、傍接円を求める。
ここでは、
\(\mathbf{p_1} = (x_1, y_1)\), \(\mathbf{p_2} = (x_2, y_2)\), \(\mathbf{p_3} = (x_3, y_3)\) \(d_1 = \|\mathbf{p_2} - \mathbf{p_3}\|\), \(d_2 = \|\mathbf{p_3} - \mathbf{p_1}\|\), \(d_3 = \|\mathbf{p_1} - \mathbf{p_2}\|\)
と定義する。
外接円 (circumcircle)
各 \(i = 1, 2, 3\) の \((x - x_i)^2 + (y - y_i)^2 = r\) から成る連立方程式を計算し、 外接円の中心 \((x, y)\) (= 外心, circumcenter) と 半径 \(r\) を求める。
1つの式から別の式を引くと以下の方程式になり \(x, y\) を求められる。
\( \displaystyle 2 \left( \begin{array}{cc} x_1 - x_2 & y_1 - y_2 \\ x_1 - x_3 & y_1 - y_3 \end{array} \right) \left( \begin{array}{c} x \\ y \end{array} \right) = \left( \begin{array}{c} (x_1^2 - x_2^2) + (y_1^2 - y_2^2) \\ (x_1^2 - x_3^2) + (y_1^2 - y_3^2) \end{array} \right)\)
内接円 (incircle)
内接円の中心 \((x, y)\) (= 内心, incenter) は以下の式で計算できる。
\(\displaystyle \left( \frac{d_1 x_1 + d_2 x_2 + d_3 x_3}{d_1 + d_2 + d_3}, \frac{d_1 y_1 + d_2 y_2 + d_3 y_3}{d_1 + d_2 + d_3} \right)\)
内接円の半径 \(r\) は \(\displaystyle \frac{2S}{d_1 + d_2 + d_3}\) で求められる。
傍接円 (excircle)
三角形に対する傍接円は3つあるが、ここでは 辺 \(\mathbf{p_2}-\mathbf{p_3}\) に対する傍接円を考える。
この傍接円の中心 \((x, y)\) (= 傍心, excenter) は以下の式で計算できる。
\(\displaystyle \left( \frac{- d_1 x_1 + d_2 x_2 + d_3 x_3}{- d_1 + d_2 + d_3}, \frac{- d_1 y_1 + d_2 y_2 + d_3 y_3}{- d_1 + d_2 + d_3} \right)\)
また、この傍接円の半径 \(r\) は \(\displaystyle \frac{2S}{- d_1 + d_2 + d_3}\) で求められる。
実装: 外接円
def circumcircle(P1, P2, P3):
x1, y1 = P1; x2, y2 = P2; x3, y3 = P3
a = 2*(x1 - x2); b = 2*(y1 - y2); p = x1**2 - x2**2 + y1**2 - y2**2
c = 2*(x1 - x3); d = 2*(y1 - y3); q = x1**2 - x3**2 + y1**2 - y3**2
det = a*d - b*c
x = d*p - b*q; y = a*q - c*p
if det < 0:
x = -x; y = -y; det = -det
x /= det; y /= det
r = ((x - x1)**2 + (y - y1)**2)**.5
return x, y, r
print("%.4f %.4f %.4f" % circumcircle((0, 0), (1, 0), (0, 1)))
# => "0.5000 0.5000 0.7071"
実装: 内接円と傍接円
def incircle(P1, P2, P3):
x1, y1 = P1; x2, y2 = P2; x3, y3 = P3
dx1 = x2 - x1; dy1 = y2 - y1
dx2 = x3 - x1; dy2 = y3 - y1
d1 = ((x3 - x2)**2 + (y3 - y2)**2)**.5
d2 = (dx2**2 + dy2**2)**.5
d3 = (dx1**2 + dy1**2)**.5
dsum = d1 + d2 + d3
r = abs(dx1 * dy2 - dx2 * dy1) / dsum
x = (x1*d1 + x2*d2 + x3*d3) / dsum
y = (y1*d1 + y2*d2 + y3*d3) / dsum
return x, y, r
def excircle(P1, P2, P3):
x1, y1 = P1; x2, y2 = P2; x3, y3 = P3
dx1 = x2 - x1; dy1 = y2 - y1
dx2 = x3 - x1; dy2 = y3 - y1
d1 = ((x3 - x2)**2 + (y3 - y2)**2)**.5
d3 = (dx1**2 + dy1**2)**.5
d2 = (dx2**2 + dy2**2)**.5
S2 = abs(dx1 * dy2 - dx2 * dy1)
dsum1 = - d1 + d2 + d3
r1 = S2 / dsum1
ex1 = (- x1*d1 + x2*d2 + x3*d3) / dsum1
ey1 = (- y1*d1 + y2*d2 + y3*d3) / dsum1
dsum2 = d1 - d2 + d3
r2 = S2 / dsum2
ex2 = (x1*d1 - x2*d2 + x3*d3) / dsum2
ey2 = (y1*d1 - y2*d2 + y3*d3) / dsum2
dsum3 = d1 + d2 - d3
r3 = S2 / dsum3
ex3 = (x1*d1 + x2*d2 - x3*d3) / dsum3
ey3 = (y1*d1 + y2*d2 - y3*d3) / dsum3
return (ex1, ey1, r1), (ex2, ey2, r2), (ex3, ey3, r3)
print("%.4f %.4f %.4f" % incircle((1, 0), (3, 3), (5, 1)))
# => "3.1472 1.5132 0.9472"
for c in excircle((1, 0), (3, 3), (5, 1)):
print("%.4f %.4f %.4f" % c)
# =>
# 5.6260 3.2600 2.0407
# 3.6726 -3.7924 4.3274
# -0.8458 2.6192 2.9887