geo2d/circle

src/geo2d/circle.h

@@ -0,0 +1,63 @@
+#pragma once
+#include "../template.h"
+#include "vector.h"
+#include "line.h" // For line intersections.
+
+// 2D circle structure; UNIT-TESTED
+//! Source: https://victorlecomte.com/cp-geo.pdf
+struct circle {
+	vec p;         // Center
+	vec::T r2 = 0; // Squared radius
+	DBP(p, r2);
+
+	// Returns -1 if point q lies outside circle,
+	// 0 if on the edge, 1 if strictly inside.
+	// Depends on vec: -, len2
+	int side(vec a) {
+		return sgn(r2 - (p-a).len2());
+	}
+
+#if FLOATING_POINT_GEOMETRY
+	// Intersect with another circle.
+	// Returns number of intersection points
+	// (greater than 2 means infinity).
+	// Arc is CCW w.r.t to `this`, CW for `a`.
+	// Depends on vec: +, -, *, len2, perp
+	int intersect(circle a, pair<vec,vec>& out) {
+		vec d = a.p - p;
+		auto d2 = d.len2();
+		if (!sgn(d2)) return !sgn(r2-a.r2) ? 3 : 0;
+		auto pd = (d2 + r2 - a.r2) / 2;
+		auto h2 = r2 - pd*pd / d2;
+		vec h, t = p + d*(pd/d2);
+		int s = sgn(h2)+1;
+		if (s > 1) h = d.perp() * sqrt(h2/d2);
+		out = {t-h, t+h};
+		return s;
+	}
+
+	// Intersect with line.
+	// Returns number of intersection points
+	int intersect(line a, pair<vec, vec>& out) {
+		auto d = a.distTo(p), h2 = r2 - d*d;
+		vec h, t = a.proj(p);
+		int s = sgn(h2)+1;
+		if (s > 1)
+			h = a.v.perp() * sqrt(h2 / a.v.len2());
+		out = {t-h, t+h};
+		return s;
+	}
+#endif
+};
+
+#if FLOATING_POINT_GEOMETRY
+// Circumcircle. Points must be non-aligned.
+// Depends on vec: +,-,*,/, cross, len2, perp
+circle circum(vec a, vec b, vec c) {
+	b = b-a; c = c-a;
+	auto s = b.cross(c);
+	assert(sgn(s));
+	vec p = a+(b*c.len2()-c*b.len2()).perp()/s/2;
+	return { p, (p-a).len2() };
+}
+#endif

tests/README.md

@@ -2,6 +2,7 @@
 
 |                                    | deterministic | fuzz            | benchmark | yosupo               |
 |------------------------------------|:-------------:|:---------------:|:---------:|:--------------------:|
+| geo2d/circle                       | partial       | partial         |           |                      |
 | geo2d/convex_hull                  |               | ✓         | ✓   |                      |
 | geo2d/halfplanes                   | ✓       | ✓         | ✓   |                      |
 | geo2d/line                         | ✓       |                 |           |                      |

tests/own/geo2d/circle.hpp

@@ -0,0 +1,116 @@
+#pragma once
+#include "../../../src/geo2d/circle.h"
+#include "common.hpp"
+
+#if FLOATING_POINT_GEOMETRY
+vec between(vec p, vec a, vec b) {
+	double x = (a-p).angle(), y = (b-p).angle();
+	if (x > y) y += M_PI*2;
+	return p + (a-p).rotate((y-x)/2);
+}
+
+void circleCircle(circle a, circle b, int expected) {
+	pair<vec, vec> out;
+	assert(a.intersect(b, out) == expected);
+
+	if (expected == 1) {
+		assert(out.x == out.y);
+		assert(a.side(out.x) == 0 && b.side(out.x) == 0);
+	} else if (expected == 2) {
+		assert(out.x != out.y);
+		assert(a.side(out.x) == 0 && b.side(out.x) == 0);
+		assert(a.side(out.y) == 0 && b.side(out.y) == 0);
+
+		// Test if arc is oriented CCW for a and CW for b.
+		assert(b.side(between(a.p, out.x, out.y)) == 1);
+		assert(a.side(between(b.p, out.y, out.x)) == 1);
+	}
+}
+
+void circleLine(circle a, line b, int expected) {
+	pair<vec, vec> out;
+	assert(a.intersect(b, out) == expected);
+
+	if (expected == 1) {
+		assert(out.x == out.y);
+		assert(a.side(out.x) == 0 && b.side(out.x) == 0);
+	} else if (expected == 2) {
+		assert(out.x != out.y);
+		assert(a.side(out.x) == 0 && b.side(out.x) == 0);
+		assert(a.side(out.y) == 0 && b.side(out.y) == 0);
+
+		// Test if points are ordered according to norm.perp()
+		assert(b.v.cross(out.x) < b.v.cross(out.y));
+	}
+}
+#endif
+
+void deterministic() {
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wfloat-equal"
+	circle def;
+	assert(def.p.x == 0 && def.p.y == 0 && def.r2 == 0);
+#pragma GCC diagnostic pop
+
+	assert(circle({4*U,2*U},25*U*U).side({1*U,7*U}) == -1);
+	assert(circle({4*U,2*U},25*U*U).side({8*U,5*U}) == 0);
+	assert(circle({4*U,2*U},25*U*U).side({7*U,5*U}) == 1);
+
+#if FLOATING_POINT_GEOMETRY
+	circleCircle({{4*U,2*U},25*U*U}, {{4*U,2*U},25*U*U}, 3);
+	circleCircle({{4*U,2*U},25*U*U}, {{4*U,2*U},16*U*U}, 0);
+	circleCircle({{4*U,2*U},25*U*U}, {{5*U,3*U},9*U*U}, 0);
+	circleCircle({{4*U,2*U},25*U*U}, {{11*U,-2*U},9*U*U}, 0);
+	circleCircle({{4*U,2*U},25*U*U}, {{11*U,-2*U},25*U*U}, 2);
+	circleCircle({{4*U,2*U},25*U*U}, {{5*U,3*U},36*U*U}, 2);
+	circleCircle({{4*U,2*U},25*U*U}, {{7*U,2*U},4*U*U}, 1);
+	circleCircle({{4*U,2*U},25*U*U}, {{7*U,2*U},4*U*U}, 1);
+	circleCircle({{4*U,2*U},25*U*U}, {{4*U,-5*U},4*U*U}, 1);
+
+	circleLine({{6*U,3*U},16*U*U}, through({1*U,2*U},{4*U,-2*U}), 0);
+	circleLine({{6*U,3*U},16*U*U}, through({2*U,7*U},{2*U,-2*U}), 1);
+	circleLine({{6*U,3*U},16*U*U}, through({3*U,7*U},{5*U,-2*U}), 2);
+	circleLine({{6*U,3*U},16*U*U}, through({5*U,-2*U},{3*U,7*U}), 2);
+#endif
+}
+
+void fuzz() {
+#if FLOATING_POINT_GEOMETRY
+	rep(i, 0, 3e6) {
+		vec p1 = randVecFromSquare(-1, 1);
+		vec p2 = randVecFromSquare(-1, 1);
+
+		int expected;
+		double d = (p1-p2).len();
+		double r2, r1 = randDouble(d/10, d*10);
+
+		if (randBool()) {
+			r2 = abs(d-r1);
+			expected = 1;
+		} else {
+			r2 = randDouble(d/10, d*10);
+			expected = (r1+r2 > d && d+r2 > r1 && d+r1 > r2 ? 2 : 0);
+		}
+
+		circle c1(p1, r1*r1), c2(p2, r2*r2);
+		if (randBool()) swap(c1, c2);
+		circleCircle(c1, c2, expected);
+	}
+
+	rep(i, 0, 3e6) {
+		vec p1 = randVecFromSquare(-1, 1);
+		vec p2 = randVecFromSquare(-1, 1);
+		vec p3 = randVecFromSquare(-1, 1);
+		auto area = (p3-p1).cross(p2-p1);
+		if (area > 1e-2) {
+			circle c = circum(p1, p2, p3);
+			assert(!sgn(c.side(p1)));
+			assert(!sgn(c.side(p2)));
+			assert(!sgn(c.side(p3)));
+		}
+	}
+#endif
+}
+
+void benchmark() {
+}

tests/own/geo2d/float/circle.cpp

@@ -0,0 +1,2 @@
+#define FLOATING_POINT_GEOMETRY 1
+#include "../circle.hpp"

tests/own/geo2d/int/circle.cpp

@@ -0,0 +1,2 @@
+#define FLOATING_POINT_GEOMETRY 0
+#include "../circle.hpp"