programming-examples/java/Data_Structures/CircleOperations.java

import java.util.*;

public class CircleOperations {

	static final double EPS = 1e-10;

	public static double fastHypot(double x, double y) {
		return Math.sqrt(x * x + y * y);
	}

	public static class Point {
		double x, y;

		public Point(double x, double y) {
			this.x = x;
			this.y = y;
		}
	}

	public static class Circle {
		double x, y, r;

		public Circle(double x, double y, double r) {
			this.x = x;
			this.y = y;
			this.r = r;
		}

		public boolean contains(Point p) {
			return fastHypot(p.x - x, p.y - y) < r + EPS;
		}
	}

	public static class Line {
		double a, b, c;

		public Line(double a, double b, double c) {
			this.a = a;
			this.b = b;
			this.c = c;
		}

		public Line(Point p1, Point p2) {
			a = +(p1.y - p2.y);
			b = -(p1.x - p2.x);
			c = p1.x * p2.y - p2.x * p1.y;
		}
	}

	// geometric solution
	public static Point[] circleLineIntersection(Circle circle, Line line) {
		double a = line.a;
		double b = line.b;
		double c = line.c + circle.x * a + circle.y * b;
		double r = circle.r;
		double aabb = a * a + b * b;
		double d = c * c / aabb - r * r;
		if (d > EPS)
			return new Point[0];
		double x0 = -a * c / aabb;
		double y0 = -b * c / aabb;
		if (d > -EPS)
			return new Point[]{new Point(x0 + circle.x, y0 + circle.y)};
		d /= -aabb;
		double k = Math.sqrt(d < 0 ? 0 : d);
		return new Point[]{
				new Point(x0 + k * b + circle.x, y0 - k * a + circle.y),
				new Point(x0 - k * b + circle.x, y0 + k * a + circle.y)};
	}

	// algebraic solution
	public static Point[] circleLineIntersection2(Circle circle, Line line) {
		return Math.abs(line.a) >= Math.abs(line.b)
				? intersection(line.a, line.b, line.c, circle.x, circle.y, circle.r, false)
				: intersection(line.b, line.a, line.c, circle.y, circle.x, circle.r, true);
	}

	static Point[] intersection(double a, double b, double c, double CX, double CY, double R, boolean swap) {
		// ax+by+c=0
		// (by+c+aCX)^2+(ay-aCY)^2=(aR)^2
		double A = a * a + b * b;
		double B = 2.0 * b * (c + a * CX) - 2.0 * a * a * CY;
		double C = (c + a * CX) * (c + a * CX) + a * a * (CY * CY - R * R);
		double d = B * B - 4 * A * C;
		if (d < -EPS)
			return new Point[0];
		d = Math.sqrt(d < 0 ? 0 : d);
		double y1 = (-B + d) / (2 * A);
		double x1 = (-c - b * y1) / a;
		double y2 = (-B - d) / (2 * A);
		double x2 = (-c - b * y2) / a;
		return swap ? d > EPS ? new Point[]{new Point(y1, x1), new Point(y2, x2)} : new Point[]{new Point(y1, x1)}
				: d > EPS ? new Point[]{new Point(x1, y1), new Point(x2, y2)} : new Point[]{new Point(x1, y1)};
	}

	public static Point[] circleCircleIntersection(Circle c1, Circle c2) {
		if (fastHypot(c1.x - c2.x, c1.y - c2.y) < EPS) {
			if (Math.abs(c1.r - c2.r) < EPS)
				return null; // infinity intersection points
			return new Point[0];
		}
		double dx = c2.x - c1.x;
		double dy = c2.y - c1.y;
		double A = -2 * dx;
		double B = -2 * dy;
		double C = dx * dx + dy * dy + c1.r * c1.r - c2.r * c2.r;
		Point[] res = circleLineIntersection(new Circle(0, 0, c1.r), new Line(A, B, C));
		for (Point point : res) {
			point.x += c1.x;
			point.y += c1.y;
		}
		return res;
	}

	public static double circleCircleIntersectionArea(Circle c1, Circle c2) {
		double r = Math.min(c1.r, c2.r);
		double R = Math.max(c1.r, c2.r);
		double d = fastHypot(c1.x - c2.x, c1.y - c2.y);
		if (d < R - r + EPS)
			return Math.PI * r * r;
		if (d > R + r - EPS)
			return 0;
		double area = r * r * Math.acos((d * d + r * r - R * R) / 2 / d / r) + R * R
				* Math.acos((d * d + R * R - r * r) / 2 / d / R) - 0.5
				* Math.sqrt((-d + r + R) * (d + r - R) * (d - r + R) * (d + r + R));
		return area;
	}

	public static Line[] tangents(Circle a, Circle b) {
		List<Line> lines = new ArrayList<>();
		for (int i = -1; i <= 1; i += 2)
			for (int j = -1; j <= 1; j += 2)
				tangents(new Point(b.x - a.x, b.y - a.y), a.r * i, b.r * j, lines);
		for (Line line : lines)
			line.c -= line.a * a.x + line.b * a.y;
		return lines.toArray(new Line[lines.size()]);
	}

	static void tangents(Point center2, double r1, double r2, List<Line> lines) {
		double r = r2 - r1;
		double z = center2.x * center2.x + center2.y * center2.y;
		double d = z - r * r;
		if (d < -EPS) return;
		d = Math.sqrt(d < 0 ? 0 : d);
		lines.add(new Line((center2.x * r + center2.y * d) / z, (center2.y * r - center2.x * d) / z, r1));
	}

	// min enclosing circle in O(n) on average
	public static Circle minEnclosingCircle(Point[] points) {
		if (points.length == 0)
			return new Circle(0, 0, 0);
		if (points.length == 1)
			return new Circle(points[0].x, points[0].y, 0);
		Collections.shuffle(Arrays.asList(points));
		Circle circle = getCircumCircle(points[0], points[1]);
		for (int i = 2; i < points.length; i++) {
			if (!circle.contains(points[i])) {
				circle = getCircumCircle(points[0], points[i]);
				for (int j = 1; j < i; j++) {
					if (!circle.contains(points[j])) {
						circle = getCircumCircle(points[j], points[i]);
						for (int k = 0; k < j; k++) {
							if (!circle.contains(points[k])) {
								circle = getCircumCircle(points[i], points[j], points[k]);
							}
						}
					}
				}
			}
		}
		return circle;
	}

	public static Circle getCircumCircle(Point a, Point b) {
		double x = (a.x + b.x) / 2.;
		double y = (a.y + b.y) / 2.;
		double r = fastHypot(a.x - x, a.y - y);
		return new Circle(x, y, r);
	}

	public static Circle getCircumCircle(Point a, Point b, Point c) {
		double Bx = b.x - a.x;
		double By = b.y - a.y;
		double Cx = c.x - a.x;
		double Cy = c.y - a.y;
		double d = 2 * (Bx * Cy - By * Cx);
		if (Math.abs(d) < EPS)
			return getCircumCircle(new Point(Math.min(a.x, Math.min(b.x, c.x)), Math.min(a.y, Math.min(b.y, c.y))),
					new Point(Math.max(a.x, Math.max(b.x, c.x)), Math.max(a.y, Math.max(b.y, c.y))));
		double z1 = Bx * Bx + By * By;
		double z2 = Cx * Cx + Cy * Cy;
		double cx = Cy * z1 - By * z2;
		double cy = Bx * z2 - Cx * z1;
		double x = cx / d;
		double y = cy / d;
		double r = fastHypot(x, y);
		return new Circle(x + a.x, y + a.y, r);
	}

	// Usage example
	public static void main(String[] args) {
		Random rnd = new Random(1);
		for (int step = 0; step < 100_000; step++) {
			int range = 10;
			int x = rnd.nextInt(range) - range / 2;
			int y = rnd.nextInt(range) - range / 2;
			int r = rnd.nextInt(range);

			int x1 = rnd.nextInt(range) - range / 2;
			int y1 = rnd.nextInt(range) - range / 2;
			int x2 = rnd.nextInt(range) - range / 2;
			int y2 = rnd.nextInt(range) - range / 2;
			if (x1 == x2 && y1 == y2)
				continue;

			Point[] p1 = circleLineIntersection(new Circle(x, y, r), new Line(new Point(x1, y1), new Point(x2, y2)));
			Point[] p2 = circleLineIntersection2(new Circle(x, y, r), new Line(new Point(x1, y1), new Point(x2, y2)));

			if (p1.length != p2.length || p1.length == 1 && !eq(p1[0], p2[0])
					|| p1.length == 2 && !(eq(p1[0], p2[0]) && eq(p1[1], p2[1]) || eq(p1[0], p2[1]) && eq(p1[1], p2[0])))
				throw new RuntimeException();
		}
	}

	static boolean eq(Point p1, Point p2) {
		return !(fastHypot(p1.x - p2.x, p1.y - p2.y) > 1e-9);
	}
}