三角形の外接円の半径と中心座標(外心)を求める

出典

Aizu Online Judge 0010
http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=0010&lang=jp

平面上の点$(x_1, y_1),(x_2, y_2), (x_3, y_3)$ を頂点とした三角形の外接円の中心座標$(p_x, p_y)$と半径rを出力するプログラムを作成して下さい。

※問題文では三角形の頂点の座標を$(x_1, y_1),(x_2, y_2), (x_3, y_3)$と定義していますが、式変形する時に訳がわからなくなるので$(a, b), (c, d), (e, f)$と読み替えます。

外心座標の計算

Wikipediaに「外心の位置」という項目がありますのでこちらが理解できる方はそのほうがたぶん早いです。私は理解できなかったので別の方法を取りました。

まず、3つの頂点からまでの外心までの距離は外接円の半径であり、全て同じであることに注目しました(下図の赤線は全て同じ長さ)。

すると、各頂点を中心とした3つの円の交点が外心座標であることに気付きます。

つまり、3つの円の方程式が並んだ3連立方程式を解けば外心座標を求められそうです。

\begin{eqnarray}
  \left\{
    \begin{array}{l}
      (p_x - a)^2 + (p_y - b)^2 = r^2　(1)\\
      (p_x - c)^2 + (p_y - d)^2 = r^2　(2)\\ 
      (p_x - e)^2 + (p_y - f)^2 = r^2　(3)\\
    \end{array}
  \right.
\end{eqnarray}

あとはゴリ押しで解くだけです。とりあえず左辺を展開します。

\begin{eqnarray}
  \left\{
    \begin{array}{l}
      p_x^2 - 2ap_x + a^2 + p_y^2 - 2bp_y + b^2 = r^2　(1')\\
      p_x^2 - 2cp_x + c^2 + p_y^2 - 2dp_y + d^2 = r^2　(2')\\
      p_x^2 - 2ep_x + e^2 + p_y^2 - 2fp_y + f^2 = r^2　(3')\\
    \end{array}
  \right.
\end{eqnarray}

$(1') - (2')$を計算すると$p_x$が求まります。

\begin{align}
  -2ap_x + 2cp_x + a^2 - c^2 - 2bp_y + 2dp_y + b^2 - d^2 &= 0 \\
  2(c - a)p_x + a^2 - c^2 + 2(d - b)p_y + b^2 -d^2   &= 0 \\
  2(c - a)p_x &= -2(d - b)p_y - a^2 - b^2 + c^2 + d^2     \\
  p_x &= \frac{2(b - d)p_y - a^2 - b^2 + c^2 + d^2}{2(c - a)}　(4)
\end{align}

同様に$(1') - (3')$を計算すると$p_x$がもう1パターン求まります。

\begin{align}
  p_x &= \frac{2(b - f)p_y - a^2 - b^2 + e^2 + f^2}{2(e - a)}　(5)
\end{align}

$(4)=(5)$なので$p_y$が求まります。

\begin{align}
  \frac{2(b - d)p_y - a^2 - b^2 + c^2 + d^2}{2(c - a)} &= \frac{2(b - f)p_y - a^2 - b^2 + e^2 + f^2}{2(e - a)} \\
  2(e - a)(b - d)p_y - (e - a)(a^2 + b^2 - c^2 - d^2) &= 2(c - a)(b - f)p_y - (c - a)(a^2 + b^2 - e^2 - f^2) \\
  2(e - a)(b - d)p_y - 2(c - a)(b - f)p_y &= (e - a)(a^2 + b^2 - c^2 - d^2) - (c - a)(a^2 + b^2 - e^2 - f^2) \\
  (2(e - a)(b - d) - 2(c - a)(b - f))p_y &= (e - a)(a^2 + b^2 - c^2 - d^2) - (c - a)(a^2 + b^2 - e^2 - f^2) \\
  p_y &= \frac{(e - a)(a^2 + b^2 - c^2 - d^2) - (c - a)(a^2 + b^2 - e^2 - f^2)}{2(e - a)(b - d) - 2(c - a)(b - f)}　(6)
\end{align}

$p_y$が求まりましたので、$(6)$を$(4)$あるいは$(5)$に代入すれば$p_x$も求まります。

まとめると以下になります。

\begin{align}
  p_y &= \frac{(e - a)(a^2 + b^2 - c^2 - d^2) - (c - a)(a^2 + b^2 - e^2 - f^2)}{2(e - a)(b - d) - 2(c - a)(b - f)} \\
  p_x &= \left\{
    \begin{array}{ll}
      \frac{2(b - d)p_y - a^2 - b^2 + c^2 + d^2}{2(c - a)} & ((c-a) \neq 0) \\
      \frac{2(b - f)p_y - a^2 - b^2 + e^2 + f^2}{2(e - a)} & ((c-a) = 0)
    \end{array}
  \right.
\end{align}

外接円の半径

外心座標$p_x$と$p_y$がわかりましたので各頂点からのユーグリッド距離を求めれば良いです。

\begin{align}
  r = \sqrt{(p_x - a)^2 + (p_y - b)^2}
\end{align}

プログラム例(Python3)

import math

n = int(input())

for i in range(n):
    points = list(map(float, input().split()))
    a = points[0]
    b = points[1]
    c = points[2]
    d = points[3]
    e = points[4]
    f = points[5]

    aa = a * a
    bb = b * b
    cc = c * c
    dd = d * d
    ee = e * e
    ff = f * f

    py = ((e - a) * (aa + bb - cc - dd) - (c - a) * (aa + bb - ee- ff)) / (2 * (e - a)*(b - d) - 2 * (c - a) * (b - f))
    
    px = (2 * (b - f) * py - aa - bb + ee + ff) / (2 * (e - a)) \
        if (c == a) else (2 * (b - d) * py - aa - bb + cc + dd) / (2 * (c - a))

    r = math.hypot(px - a, py - b)

    print("%.3f %.3f %.3f" % (px, py, r))

蛇足

なんだかゴリ押し過ぎる。$p_x$を求めるにも条件分岐が必要だし美しくないですね・・・他にもっときれいな解法がありそう。
まあ解けたのでよしとするか。