三角形の外接円の中心を求めるコード

ただなんとなく学校のJavaの課題のコードがあまりにも計算機パワーに頼りがちなコードだったので、これを 外接円の中心は２つの垂直二等分線の交点 という基礎的な事実に基づいてコードを書いてみた。

ちなみに a2*a6-a3*a5 が0に近づいたときなどのエラー処理はおこなってない。

あと、 と で、とが直交するとき、 になる理屈がちょっと分からなくて混乱した。まあ、を使えば簡単に分かったけど。

circle.c （私が書いたやつ）

#include <stdio.h>

typedef struct __Point {
	double x;
	double y;
} Point;

void CalcPoint(Point p, Point q, Point r, Point *ret)
{
	double a1 = q.y + r.y;
	double a4 = q.x + r.x;
	double a7 = p.x + q.x;
	double a8 = p.y + q.y;
	double a2 = q.x - r.x;
	double a3 = q.y - r.y;
	double a5 = p.x - q.x;
	double a6 = p.y - q.y;
	
	ret->x = (a1*a3*a6+a2*a4*a6-a3*a5*a7-a3*a6*a8)/(a2*a6-a3*a5)/2.0;
	ret->y = -a5/a6 * ret->x + a8/2.0 + (a5*a7)/(a6*2.0);
	
	return;
}

int main(int argc, char *argv[])
{
	Point p;
	Point q;
	Point r;
	
	Point ans;
	
	/* 点(15, 1) */
	p.x = 15, p.y = 1;
	
	/* 点(1, 5) */
	q.x = 1, q.y = 5;
	
	/* 点(4, 0) */
	r.x = 4, r.y = 0;
	
	CalcPoint(p, q, r, &ans);
	
	printf("x = %.15lf, y = %.15lf\n", ans.x, ans.y);

	return 0;
}

また、学校で使われたコードを載せてみようと思う。Javaなので関係するメソッドだけを抜粋する。それにしても、あんまりだよなあ、と感じる。

Point.java （一部分）

	/**
	 * 自分自身と引数で指定された2つの点の合計3点から、
	 * ほぼ同一の距離にある点を近似的に計算して求める
	 * 
	 * calculate approximate center point of this, p1 and p2
	 * 
	 * @param p1
	 *            destination1
	 * @param p2
	 *            destination2
	 * @return approximate center point
	 */
	public Point center(Point p1, Point p2) {
		// define candidate point
		Point p = center(p1).center(p2);
		Point[] pp = new Point[] { this, p1, p2 };
		double[] dist = new double[pp.length];

		int count = 0;
		double maxdiff;
		do {
			// calculate distances from candidate to 3 points
			double avg = 0;
			for (int i = 0; i < pp.length; i++) {
				dist[i] = p.distanceTo(pp[i]);
				avg += dist[i];
			}
			avg /= pp.length;

			// edge up
			double dx = 0.0, dy = 0.0, dz = 0.0;
			maxdiff = 0.0;
			for (int i = 0; i < pp.length; i++) {
				dist[i] /= avg;
				if (maxdiff < Math.abs(dist[i] - 1.0)) {
					maxdiff = Math.abs(dist[i] - 1.0);
				}
				dx += (Math.sqrt(dist[i]) - 1.0) * (pp[i].x - p.x);
				dy += (Math.sqrt(dist[i]) - 1.0) * (pp[i].y - p.y);
				dz += (Math.sqrt(dist[i]) - 1.0) * (pp[i].z - p.z);
			}
			p.x += dx;
			p.y += dy;
			p.z += dz;

			count++;
		} while (maxdiff > 1.0E-15 && count < 3000);

		return p;
	}